Tryfon Antonakakis

Fadi Baïda

Abderrahmane Belkhir

Kirill Cherednichenko

Shane Cooper

Richard Craster

Guillaume Demesy

John Desanto

Gérard Granet

Editor Evgeny Popov

Boris Gralak
email: boris.gralak@fresnel.fr

Sébastien Guenneau

Daniel Maystre

André Nicolet

Brian Stout

Fréderic Zolla

Benjamin Vial

GratinGs: theory and numeric applications

Introduction

Exact modal method (EMM) has been proposed to take advantage of geometry of lamellar gratings. These gratings are made of rectangular rods periodically spaced which can be considered locally as periodic multilayered stacks (see figure 10.1). This simple geometry makes it possible to expand the electromagnetic field on the basis of "exact modes", and to obtain an exact representation of the permittivity. In this particular case, EMM can be more efficient than similar methods based on Fourier expansion (coupled-wave method [START_REF] Moharam | Rigorous coupled-waves analysis of metallic surfacerelief grating[END_REF] or Fourier modal method [START_REF] Li | New formulation of the Fourier modal method for crossed surface-relief gratings[END_REF][START_REF] Li | Justification of matrix truncation in the modal methods of diffraction gratings[END_REF]) which may lead to poor convergence due to the discontinuous nature of both electromagnetic field and permittivity. This advantage of EMM becomes more important when the permittivity contrast is high, e.g. for metallic lamellar gratings. . . . . . .

.1: A lamellar grating made of a single lamellar layer on a substrate. The region corresponding to the lamellar layer, between the planes x 3 = 0 and x 3 = h, can be considered as the multilayered stack on the left.

Exact modal method has been introduced in 1981 in order to solve Maxwell's equations in presence of lamellar gratings made of dielectrics [START_REF] Botten | The dielectric lamellar diffraction grating[END_REF] and metals [START_REF] Botten | The finitely conducting lamellar diffraction grating[END_REF][START_REF] Botten | Highly conducting lamellar diffraction grating[END_REF][START_REF] Tayeb | On the numerical study of deep conducting lamellar diffraction grating[END_REF]. Since these pioneering works, a major contribution to this method is certainly its rigorous extension to conical mountings [START_REF] Li | A modal analysis of lamellar diffraction gratings in conical mountings[END_REF], on which is based an EMM for three-dimensional woodpile structures [START_REF] Gralak | Theoretical study of photonic band gaps in woodpile crystals[END_REF]. Another major development is the introduction of perfectly matched layers in order to model aperiodic systems met in integrated optics [START_REF] Silberstein | Use of grating theories in integrated optics[END_REF] (information can be found on the website of CAMFR).

In this chapter, a rigorous formulation of the exact modal method for lamellar structures is presented. In section 10.3, a special attention is paid to the continuation of the electromagnetic field inside a lamellar layer. In combination with the boundary conditions, this continuation provides a large class of solutions of Maxwell's equations in presence of lamellar gratings. In section 10.4, it is shown that, in each lamellar layer, there is a decoupling of the vector field equations into two independent scalar equations, which correspond to the ones of a multilayered stack (see figure 10.1). Numerical stacking algorithms are presented in section 10.5 and a numerical illustration of the EMM efficiency is proposed in section 10.6. Finally, the techniques used for the calculation of the exact modes and the associated exact eigenvalues are reported in the appendix (section 10.7).

Note that the important extensions to woodpile structures [START_REF] Gralak | Theoretical study of photonic band gaps in woodpile crystals[END_REF] and to lamellar gratings including infinitely conducting metal [START_REF] Gralak | Solutions of Maxwell's equations in presence of lamellar gratings including infinitely conducting metal[END_REF] are not considered in this book chapter. These cases will be however included in the next version.

Notations

Throughout this chapter an orthonormal basis (e e e 1 , e e e 2 , e e e 3 ) is used: every vector x x x in R 3 is described by its three components x 1 , x 2 and x 3 . It is shown how to obtain in the presence of a stack of lamellar layers, a large class of solutions E E E of the Helmholtz equation

ω 2 -ε -1 ∇ ∇ ∇ × µ -1 ∇ ∇ ∇ × E E E = 0 0 0 , (10.1)
where ε is the permittivity, µ is the permeability, ω is the frequency and ∇ ∇ ∇× is the curl operator.

All the media considered in this chapter are isotropic, and thus the permittivity and permeability reduce to scalar functions. The considered structure is independent of the variable x 2 , and x 1periodic with spatial period d d d = de e e 1 :

ε(x x x + d d d) = ε(x x x) = ε(x 1 , x 3 ) , µ(x x x + d d d) = µ(x x x) = µ(x 1 , x 3 ) , x x x ∈ R 3 . (10.
2)

The unit cell associated with this grating is [0, d] and the one-dimensional lattice is nd n ∈ Z . Then, a lamellar grating is a stack in the direction x 3 of lamellar layers where ε and µ are both functions of the single variable x 1 (figure 10.2). In practice, each lamellar layer is made of infinite parallel rods with rectangular cross section (figure 10.2): the functions ε and µ are piecewise constant of the solely variable x 1 .

x 3

x 1

x 2 d In order to obtain a set of first order differential equations from (10.1) a second field is defined:

H H H = (ω µ) -1 ∇ ∇ ∇ × E E E . (10.3)
Note that this quantity differs from the usual "harmonic H H H field" by the complex number i.

Solutions E E E, H H H are investigated in the space of fields whose restrictions in every horizontal plane (normal to e e e 3 ) are square integrable:

R 2 F F F(x 1 , x 2 , x 3 ) 2 dx 1 dx 2 < ∞ , x 3 ∈ R , (10.4)
where

F F F = E E E, H H H.
The first consequence of (10.4) is the possibility to perform a decomposition of the problem to take advantage of the spatial invariances of the system: a Fourier decomposition with respect to the variable x 2 ,

F F F -→ F F F(x 1 , k 2 , x 3 ) = 1 2π R exp[-ik 2 x 2 ]F F F(x 1 , x 2 , x 3 ) dx 2 , (10.5)
and a Floquet-Bloch decomposition with respect to the variable x 1 ,

F F F -→ F F F(k 1 , x 1 , k 2 , x 3 ) 1 2π ∑ n∈Z exp[-ik 1 nd] F F F(x 1 + pd, k 2 , x 3 ) , (10.6)
where k 1 is the Bloch wave vector in the first Brillouin zone

[-π/d, π/d]. Thus solutions Ẽ E E, H H H satisfy [-π/d,π/d] F F F(k 1 , x 1 , k 2 , x 3 ) 2 dx 1 < ∞ , x 1 , k 2 , x 3 ∈ R , (10.7)
with the partial Bloch boundary condition

F F F(k 1 , x 1 + d, k 2 , x 3 ) = exp[ik 1 d] F F F(k 1 , x 1 , k 2 , x 3 ) , (10.8)
where

k 1 is fixed in [-π/d, π/d].
The second consequence of (10.4) [or (10.7)] is that the restrictions to every horizontal plane of ∇ ∇ ∇ × E E E and ∇ ∇ ∇ × H H H are also locally square integrable [from (10.1,10.3)]. Then, for all i, j = 1, 2, 3 and i = j, E i and H i are continuous functions of the variable x j . In particular, the tangential components E 1 , E 2 , H 1 and H 2 of E E E and H H H are continuous functions of the variable x 3 . It follows that it is possible to solve Maxwell's equations in a stack of layers by the following two steps: the first step consists in solving Maxwell's equations in each layer independently and then the second step consists in connecting each independent solution using the continuity of

E 1 , E 2 , H 1 and H 2 .
With the definition (10.3), equation (10.1) is equivalent to the set of first order equations

E E E = (ωε) -1 ∇ ∇ ∇ × H H H , H H H = (ω µ) -1 ∇ ∇ ∇ × E E E . (10.9)
Let the 2 × 2 matrix σ and the two-components vector F j defined by

σ = ω 0 µ ε 0 , F j = Ẽ j Hj , j = 1, 2, 3 . (10.10)
Then, the first order equations (10.9) can be developed as

F 1 = σ -1 ∂ 2 F 3 -∂ 3 F 2 , F 2 = σ -1 ∂ 3 F 1 -∂ 1 F 3 , F 3 = σ -1 ∂ 1 F 2 -∂ 2 F 1 , (10.11)
where ∂ j is the partial derivative with respect to the variable x j ( j = 1, 2, 3). This last set of equations is exactly the same as (10.9) and, with some abuse of notations, it can written in the compact way

F F F = σ -1 ∇ ∇ ∇ × F F F, with F F F = (F 1 , F 2 , F 3 ).

Continuation of the electromagnetic field

In this section two different formulations are presented to solve the equation F F F = σ -1 ∇ ∇ ∇ × F F F in a lamellar layer located between the planes x 3 = 0 and x 3 = h. In practice, this solution is expressed as a relationship between F F F(0) and F F F(h). Note that the formulations presented in this section remain valid in the general case of cross gratings with two-dimensional periodicity [START_REF] Li | New formulation of the Fourier modal method for crossed surface-relief gratings[END_REF].

Direct formulation: the transfer matrix

The starting point is the set of equations (10.11). Eliminating the components F 3 , one obtains

∂ 3 F = iM F , F = F 1 F 2 , M = -i -∂ 1 σ -1 ∂ 2 σ + ∂ 1 σ -1 ∂ 1 -σ -∂ 2 σ -1 ∂ 2 ∂ 2 σ -1 ∂ 1 . (10.12)
For a lamellar layer located between the planes x 3 = 0 and x 3 = h (see figure 10.1), the functions ε and µ are x 3 -independent for x 3 in [0, h]: then the matrix σ and the operator-valued matrix M are also x 3 -independant for x 3 in [0, h]. Let L be the x 3 -independant operator valued matrix which coincides with M in this single layer:

L = M(x 3 ) , x 3 ∈ [0, h] . (10.13)
In a first step, it is assumed that (see the end of section 10.4 for a justification) the matrix L has a diagonal form and can be written as

L = V λ V -1 , (10.14)
where the matrix V contains the eigenvectors V ±,n of L, and λ is the diagonal matrix made of the associated eigenvalues λ ±,n :

LV ±,n = λ ±,n V ±,n . (10.15)
The sets of eigenvectors and eigenvalues are split into two parts according to the sign of the imaginary part of λ ±,n : Imλ +,n > 0 and Imλ -,n < 0. Let λ + (respectively λ -) be the diagonal matrices containing the eigenvalues λ +,n > 0 of L (respectively λ -,n > 0): then .16) This last condition on the imaginary part of eigenvalues is always realized if there is some absorption, i.e. Imσ > 0, or if a small positive imaginary part is added to the frequency ω (in the later case the limit ω = lim η↓0 (ω + iη) is considered [START_REF] Tip | Band structure of absorptive photonic crystals[END_REF][START_REF] Gralak | Transfer matrix method for point sources radiating in classes of negative refractive index materials with 2n-fold antisymmetry[END_REF]). The combination of (10.12) and (10.13) leads to the equation

λ = λ + 0 0 λ - , Imλ + > 0 , Imλ -< 0 . ( 10
∂ 3 F(x 3 ) = iL F(x 3 ) , x 3 ∈ [0, h] , (10.17)
where the dependence on other variables has been omitted. Since L is x 3 -independent, the "formal" solution of this equation is just .18) This solution is denominated by "formal" since, at this stage, it is still necessary to check if it exists. Using the diagonal form (10.14) of the operator L, the expression (10.18) becomes

F(x 3 ) = exp[iLx 3 ]F(0) , x 3 ∈ [0, h] . ( 10
F(x 3 ) = V exp[iλ x 3 ]V -1 F(0) . (10.19)
Actually, the diagonal matrix exp[iλ x 3 ] is made of the two parts exp[iλ ± x 3 ] which have different behaviour. From (10.16), the part exp[iλ + x 3 ] is bounded by exp[-Imλ + x 3 ] < 1. On the contrary, the part exp[iλ -x 3 ] is not bounded and, in general, the corresponding coefficients are growing towards infinity like exponential functions. Consequently, the transfer matrix T (x 3 ) defined by

F(x 3 ) = T (x 3 )F(0) , T (x 3 ) = V exp[iλ x 3 ]V -1 , (10.20)
has "infinite" coefficients and thus expressions like (10.18), (10.19) and (10.20) have to be considered as purely "formal" and have to be handled cautiously. Numerically, the transfer matrix is truncated and, because its coefficients tend to infinity like exponential functions, it presents numerical instabilities which makes it difficult to use it. A numerical solution has been found to solve this problem with the definition of the Sand R-algorithms [START_REF] Li | Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings[END_REF] (see section 10.5 for the numerical solution in the present case). Finally, notice that the transfer matrix is occasionally derived from the matrix L 2 instead of L. Indeed, the assumption (10.14) on the diagonal form of L might be too strong and not rigorously true. According to notations (10.14), we denote by V and λ 2 the matrices containing the eigenvectors and eigenvalues of L 2 :

L 2 = V λ 2 V -1 . (10.21)
In that case, it used that equations (10.12) and (10.13) imply

∂ 2 3 F = -L 2 F , x 3 ∈ [0, h] . (10.22)
The combination of the two last equations leads to

F(x 3 ) = V cos[λ x 3 ]V -1 F(0) +V λ -1 sin[λ x 3 V -1 (∂ 3 F)(0) . (10.23)
Replacing (∂ 3 F)(0) by iLF(0), one obtains for the transfer matrix the following "formal" expression

T (x 3 ) = V cos[λ x 3 ]V -1 + iV λ -1 sin[λ x 3 ]V -1 L . (10.24)
This equation is not "formally" equivalent to the first expression (10.20) derived from L. This equivalence requires the assumption (10.14) to become true, so that L can be replaced by V λV -1 above (and next the formal identity cos[λ

x 3 ] + i sin[λ x 3 ] = exp[iλ x 3 ]
has to be used).

Rigorous derivation of the continuation procedure

A rigorous formulation is based on the use of the Fourier transform with respect to the variable x 3 defined by

F [F](k 3 ) = 1 2π R exp[-ik 3 x 3 ]F(x 3 ) dx 3 . (10.25)
The function F is then deduced from its Fourier transform F [F] by

F(x 3 ) = R exp[ik 3 x 3 ]F [F](k 3 ) dk 3 . (10.26)
It is not suitable to perform directly the Fourier transform of the equation (10.12) since the matrix σ (and then M) is not independent of x 3 in R. However, if equation (10.12) is multiplied by the characteristic function Ψ of the lamellar layer [so Ψ(x 3 ) = 1 for x 3 in [0, h] and vanishes otherwise], then

Ψ(x 3 )∂ 3 F(x 3 ) = Ψ(x 3 )iM F(x 3 ) = Ψ(x 3 )iL F(x 3 ) = iL Ψ(x 3 )F(x 3 ) . (10.27)
After this multiplication, a partial differential equation with the x 3 -independent matrix L is obtained. The Fourier transform (10.25) of Ψ∂ 3 F is

F [Ψ∂ 3 F](k 3 ) = 1 2π R exp[-ik 3 x 3 ]Ψ(x 3 )∂ 3 F(x 3 ) dx 3 = 1 2π h 0 exp[-ik 3 x 3 ]∂ 3 F(x 3 ) dx 3 = 1 2π exp[-ik 3 h]F(h) -F(0) + ik 3 F [ΨF](k 3 ) , (10.28)
where the last line comes from an integration by parts. After this Fourier transform, equation (10.27) becomes

1 2π exp[-ik 3 h]F(h) -F(0) + ik 3 F [ΨF](k 3 ) = iLF [ΨF](k 3 ) (10.29) or k 3 -L F [ΨF](k 3 ) = 1 2iπ F(0) -exp[-ik 3 h]F(h) . (10.30)
The operator k 3 -L is always invertible if there is some absorption, i.e. Imσ > 0, or if the limit ω = lim η↓0 (ω + iη) is considered (see [START_REF] Tip | Band structure of absorptive photonic crystals[END_REF][START_REF] Gralak | Transfer matrix method for point sources radiating in classes of negative refractive index materials with 2n-fold antisymmetry[END_REF], this is equivalent to the property (10.16) on the eigenvalues λ since k 3 is purely real). Hence it is possible to write

F [ΨF](k 3 ) = 1 2iπ 1 k 3 -L F(0) -exp[-ik 3 h]F(h) , (10.31)
The final step is to apply the inverse Fourier transform (10.26): for x 3 in [0, h],

Ψ(x 3 )F(x 3 ) = 1 2iπ R exp[ik 3 x 3 ] 1 k 3 -L dk 3 F(0) - 1 2iπ R exp[ik 3 (x 3 -h)] 1 k 3 -L dk 3 F(h) . (10.32)
Again, it is assumed that the operator L can be written L = V λV -1 (10.14). Replacing the matrix L by its diagonal form, the last expression becomes

Ψ(x 3 )F(x 3 ) = 1 2iπ V R exp[ik 3 x 3 ] 1 k 3 -λ dk 3 V -1 F(0) - 1 2iπ V R exp[ik 3 (x 3 -h)] 1 k 3 -λ dk 3 V -1 F(h) .
(10.33)

The integrations above are performed by adding to the real axis of k 3 a semi-circle with infinite radius (in the complex plane of k 3 ) on which the integrals vanish. For the first term with F(0), the complex number k 3 must have positive imaginary part (x 3 is positive), so that the real axis is closed by a semi-circle in the upper half plane (see the red path on figure 10.3). In this case, the solely eigenvalues contained in λ + generate contributions in the integral. For the second term with F(h), the complex number k 3 must have negative imaginary part (x 3h is negative), so that the real axis is closed by a semi-circle in the lower half plane (see the blue path on figure 10.3). Here, the integral is given by the eigenvalues contained in λ -. Let P ± be the projectors Re(k 3 ) upon the spaces corresponding respectively to eigenvalues λ ± :

Im(k 3 ) exp[ik 3 x 3 ] → 0 exp[ik 3 (x 3 -h)] → 0 λ + λ -
P + λ = λ + 0 0 0 , P -λ = 0 0 0 λ - . (10.34)
Then, after the integration over k 3 , expression (10.33) yields

Ψ(x 3 )F(x 3 ) = V P + exp[iλ + x 3 ]V -1 F(0) +V P -exp[iλ -(x 3 -h) V -1 F(h) . (10.35)
This expression is always well defined since the integration in the complex plane of k 3 imposes that all the complex exponential functions decrease:

exp[iλ + x 3 ] ≤ 1 , exp[iλ -(x 3 -h)] ≤ 1 . (10.36)
Considering the rigorous expression (10.35) at x 3 = 0 and x 3 = h and using that P -+ P + is the identity, one obtains

V P -V -1 F(0) = V P -exp[-iλ -h]V -1 F(h) , V P + V -1 F(h) = V P + exp[iλ + h]V -1 F(0) . (10.37)
These two relationships provides a rigorous way to deduce F(0) from F(h) and conversely. As in the previous section, the continuation procedure is also derived from the diagonal form (10.21) of L 2 . Equation (10.22) is multiplied by Ψ(x 3 ) to provide an expression similar to (10.27)

Ψ(x 3 )∂ 3 F(x 3 ) = -L 2 Ψ(x 3 )F(x 3 ) . (10.38)
Then, the Fourier transform (10.25) is applied to this equation. Using that

F [Ψ∂ 2 3 F](k 3 ) = 1 2π R exp[-ik 3 x 3 ]Ψ(x 3 )∂ 2 3 F(x 3 ) dx 3 = -k 2 3 F [ΨF](k 3 ) + ik 3 1 2π exp[-ik 3 h]F(h) -F(0) + 1 2π exp[-ik 3 h](∂ 3 F)(h) -(∂ 3 F)(0) , (10.39)
and replacing (∂ 3 F)(x 3 ) by iLF(x 3 ), equation ( 10.38) implies

F [ΨF](k 3 ) = 1 2iπ 1 k 2 3 -L 2 k 3 F(0) -exp[-ik 3 h]F(h) = 1 2iπ 1 k 2 3 -L 2 LF(0) -exp[-ik 3 h]LF(h) .
(10.40)

Next, the inverse Fourier transform (10.26) is performed for x 3 in [0, h] and the diagonal form (10.21) is used:

Ψ(x 3 )F(x 3 ) = 1 2iπ V R k 3 k 2 3 -λ 2 exp[ik 3 x 3 ]V -1 F(0) -exp[ik 3 (x 3 -h)]V -1 F(h) dx 3 + 1 2iπ V R 1 k 2 3 -λ 2 exp[ik 3 x 3 ]V -1 LF(0) -exp[ik 3 (x 3 -h)]V -1 LF(h) dx 3 .
(10.41) Again, the integrations above are calculated by adding to the real axis of k 3 a semi-circle with infinite radius (in the complex plane of k 3 ) on which the integrals vanish. Without loss of generality, it is considered that the square root of the eigenvalues in λ 2 have non-zero imaginary part: let √ λ 2 be the square root of λ 2 with positive imaginary part. For the terms with F(0), the real axis is closed by a semi-circle in the upper half plane, and the eigenvalues with positive imaginary part √ λ 2 lead to contributions in the integrals. For the terms with F(h), the real axis is closed by a semi-circle in the lower half plane, and the integrals are given by the eigenvalues with negative imaginary part, i.e. -√ λ 2 . Calculations of integrals over k 3 lead to

Ψ(x 3 )F(x 3 ) = V 1 2 exp[i √ λ 2 x 3 ]V -1 F(0) +V 1 2 exp[-i √ λ 2 (x 3 -h)]V -1 F(h) + V 1 2 √ λ 2 exp[i √ λ 2 x 3 ]V -1 LF(0) -V 1 2 √ λ 2 exp[-i √ λ 2 (x 3 -h)]V -1 LF(h) .
(10.42) This equation is evaluated at x 3 = 0 and x 3 = h:

F(0) = V exp[i √ λ 2 h]V -1 F(h) +V 1 √ λ 2 V -1 LF(0) -V 1 √ λ 2 exp[i √ λ 2 h]V -1 LF(h) , F(h) = V exp[i √ λ 2 h]V -1 F(0) +V 1 √ λ 2 exp[i √ λ 2 h]V -1 LF(0) -V 1 √ λ 2
V -1 LF(h) .

(10.43) Thanks to the technique based on the Fourier transform, all the exponential functions in these expressions must be well-defined. Indeed, the imaginary part of √ λ 2 is positive and all the exponential functions decrease. Equations (10.43) will be used to construct a stable numerical algorithm to stack several lamellar layers.

Exact eigenmodes and eigenvalues method

The different solutions (10.20), (10.24, (10.37) and (10.43), established in the previous section, are provided from the knowledge of the sets of eigenmodes and eigenvalues of the operator L (or L 2 ). In this section, it is shown how these eigenmodes V ±,n and eigenvalues λ ±,n of L can be exactly determined in a lamellar layer located between the planes x = 0 and x 3 = h. The starting point is equation (10.17):

∂ 3 F = iLF , L = M(x 3 ) , x 3 ∈ [0, h]. (10.44)
Let ε 1 , µ 1 and σ 1 be the functions which coincide with respectively ε, µ and σ in the considered lamellar layer for x 3 in [0, h]. In the considered lamellar layer, they are functions of the solely variable x 1 (see figure 10.2):

ε 1 (x 1 ) = ε(x 1 , x 3 ) , µ 1 (x 1 ) = µ(x 1 , x 3 ) , σ 1 (x 1 ) = σ (x 1 , x 3 ) , x 3 ∈ [0, h]. (10.45)
According to expression (10.12), the operator L is now

L = -i -∂ 1 σ -1 1 ∂ 2 σ 1 + ∂ 1 σ -1 1 ∂ 1 -σ 1 -∂ 2 σ -1 1 ∂ 2 ∂ 2 σ -1 1 ∂ 1 , (10.46)
and its square is

L 2 = - -σ 2 1 -∂ 1 σ -1 1 ∂ 1 σ 1 -σ 1 ∂ 2 σ -1 1 ∂ 2 σ 1 ∂ 2 σ -1 1 ∂ 1 -∂ 1 σ -1 1 ∂ 2 σ 1 σ 1 ∂ 1 σ -1 1 ∂ 2 -∂ 2 σ -1 1 ∂ 1 σ 1 -σ 2 1 -∂ 2 σ -1 1 ∂ 2 σ 1 -σ 1 ∂ 1 σ -1 1 ∂ 1 . (10.47)
Since the matrix σ 1 is x 2 -independent, the equality

σ 1 ∂ 2 σ -1 1 = ∂ 2 = σ -1 1 ∂ 2 σ 1

holds, and the expression above becomes

L 2 = σ 2 1 + ∂ 2 2 + ∂ 1 σ -1 1 ∂ 1 σ 1 0 ∂ 2 σ -1 1 ∂ 1 σ 1 -σ 1 ∂ 1 σ -1 1 ∂ 2 σ 2 1 + ∂ 2 2 + σ 1 ∂ 1 σ -1 1 ∂ 1 . (10.48)
This expression shows that the components F 1 can be decoupled from the components F 2 in the lamellar layer. Indeed, it implies

∂ 2 3 F 1 = -KF 1 , K = σ 2 1 + ∂ 2 2 + ∂ 1 σ -1 1 ∂ 1 σ 1 . (10.49)
Moreover, each component of F 1 , i.e. E 1 and H 1 , can be also decoupled since the operator K is diagonal:

K = K ε 1 0 0 K µ 1 , (10.50)
where

K ε 1 = ω 2 ε 1 µ 1 + ∂ 2 2 + ∂ 1 ε -1 1 ∂ 1 ε 1 , K µ 1 = ω 2 ε 1 µ 1 + ∂ 2 2 + ∂ 1 µ -1 1 ∂ 1 µ 1 .
(10.51)

Here, it is important to notice that the two operators K ε 1 and K µ 1 correspond to the ones of a onedimensional multilayered stack for respectively pand s-polarization. This makes it possible to calculate the exact eigenmodes and eigenvalues of K ε 1 and K µ 1 (see appendix) and thus the ones of K. Thus the continuation procedure presented in section 10.3.2 can be applied to equation (10.49). It provides relationships between the fields F 1 (0), F 1 (h) and their derivative with respect to x 3 , i.e.

[∂ 3 F 1 ](0) and [∂ 3 F 1 ](h) (it is recalled that ∂ 3 F = iLF in section 10.3.2).
To complete the derivation of the method, it is necessary to express the component F 2 of the field from F 1 and ∂ 3 F 1 . A starting relationship is obtained from (10.44) and (10.47):

∂ 3 F 1 = -∂ 1 σ -1 1 ∂ 2 F 1 + σ 1 + ∂ 1 σ -1 1 ∂ 1 F 2 . (10.52)
This equation is equivalent to

1 + ∂ 1 σ -1 1 ∂ 1 σ -1 1 σ 1 F 2 = ∂ 3 F 1 + ∂ 1 σ -1 1 ∂ 2 F 1 . (10.53)
Here, it is remarked that the operator

[1 + ∂ 1 σ -1 1 ∂ 1 σ -1 1 ] is invertible since it equals [K -∂ 2 2 ]σ -2 1
where σ 1 is invertible as well as [K -∂ 2 2 ] (eigenvalues of K have non-zero imaginary part). Moreover, the operator

[1 + ∂ 1 σ -1 1 ∂ 1 σ -1 1 ] commutes with ∂ 3 , ∂ 2 and ∂ 1 σ -1 1 : hence σ 1 F 2 = 1 1 + ∂ 1 σ -1 1 ∂ 1 σ -1 1 ∂ 3 F 1 + 1 1 + ∂ 1 σ -1 1 ∂ 1 σ -1 1 ∂ 1 σ -1 1 ∂ 2 F 1 = 1 1 + ∂ 1 σ -1 1 ∂ 1 σ -1 1 ∂ 3 F 1 + ∂ 2 ∂ 1 σ -1 1 1 1 + ∂ 1 σ -1 1 ∂ 1 σ -1 1 F 1 = σ 2 1 1 σ 2 1 + ∂ 1 σ -1 1 ∂ 1 σ 1 ∂ 3 F 1 + ∂ 2 ∂ 1 σ 1 1 σ 2 1 + ∂ 1 σ -1 1 ∂ 1 σ 1 F 1 .
(10.54)

After the Fourier decomposition with respect to the variable x 2 , the operator ∂ 2 becomes ik 2 .

Then, the inverse operator above is expressed using the diagonal form K = V λ 2 V -1 , and the component F 2 becomes

F 2 = σ 1 V 1 λ 2 + k 2 2 V -1 ∂ 3 F 1 + ik 2 σ -1 1 ∂ 1 σ 1 V 1 λ 2 + k 2 2 V -1 F 1 . (10.55)
Finally, it is stressed that the coefficients of the operators σ 1 V and σ -1 1 ∂ 1 σ 1 V above can be calculated exactly from the knowledge of the exact eigenmodes. Indeed, the technique presented in appendix shows that the determination of the eigenmodes φ n lies on the calculations of functions σ 1 φ n and σ -1

1 ∂ 1 σ 1 φ n . The final expression of F 2 in terms of F 1 and ∂ 3 F 1 is then

F 2 = U 1 λ 2 + k 2 2 V -1 ∂ 3 F 1 + ik 2 W 1 λ 2 + k 2 2 V -1 F 1 , U = σ 1 V , W = σ -1 1 ∂ 1 σ 1 V . ( 10
.56) To conclude this section, it has been shown that the field components F 1 , ∂ 3 F 1 and F 2 can be expressed exactly in a lamellar layer: F 1 and ∂ 3 F 1 are provided by equations (10.42) and (10.43) (where F → F 1 , LF → -i∂ 3 F 1 , and K = V λ 2 V -1 ); F 2 is given by equation (10.56). These expressions are based on the knowledge of the eigenvectors V of K, the eigenvalues λ 2 of K and the operators U = σ 1 V and W = σ -1 1 ∂ 1 σ 1 V . All these quantities can be calculated exactly, as shown in appendix.

Finally, it is stressed that the propagation constants inside the lamellar layer are the square roots of the eigenvalues of the operator K, i.e. ± √ λ 2 . If there is some absorption, i.e. Imσ > 0, or if a small positive imaginary part is added to the frequency ω (the limit ω = lim η↓0 (ω + iη) is considered [START_REF] Tip | Band structure of absorptive photonic crystals[END_REF][START_REF] Gralak | Transfer matrix method for point sources radiating in classes of negative refractive index materials with 2n-fold antisymmetry[END_REF]), the eigenvalues λ 2 of K cannot be purely real. It follows that the propagation constants ± √ λ 2 have non zero imaginary parts and, moreover, the half of them (+ √ λ 2 ) have strictly positive imaginary part and the second half of them (-√ λ 2 ) have strictly negative imaginary part. In this case, the assumption (10.16) on the eigenvalues of L is well justified.

Numerical algorithm

A general solution for numerical algorithm has been proposed by L. Li in the case of modal methods of gratings [START_REF] Li | Justification of matrix truncation in the modal methods of diffraction gratings[END_REF]. This solution is based on the definition of S or R matrices which are well-conditioned.

R matrix for a single lamellar layer

In this section, the expression of a R matrix associated with a lamellar layer is established: it is defined by the relationship

F 1 (0) F 1 (h) = R F 2 (0) F 2 (h) . (10.57)
First, equation ( 10.43) is used to provide an expression of the solution of (10.49):

F 1 (0) -V exp[i √ λ 2 h]V -1 F 1 (h) = -iV 1 √ λ 2 V -1 ∂ 3 F 1 (0) + iV 1 √ λ 2 exp[i √ λ 2 h]V -1 ∂ 3 F 1 (h) , F 1 (h) -V exp[i √ λ 2 h]V -1 F 1 (0) = -iV 1 √ λ 2 exp[i √ λ 2 h]V -1 ∂ 3 F 1 (0) + iV 1 √ λ 2 V -1 ∂ 3 F 1 (h) .
(10.58) This set of equations is written using 2 × 2 matrices:

A F 1 (0) F 1 (h) = B ∂ 3 F 1 (0) ∂ 3 F 1 (h) , (10.59)
where

A = 1 -V exp[i √ λ 2 h]V -1 -V exp[i √ λ 2 h]V -1 1 (10.60) and B =     -iV 1 √ λ 2 V -1 iV 1 √ λ 2 exp[i √ λ 2 h]V -1 -iV 1 √ λ 2 exp[i √ λ 2 h]V -1 iV 1 √ λ 2 V -1     .
(10.61)

Next, from (10.56), the field ∂ 3 F 1 is related to F 1 and F 2 from

∂ 3 F 1 = -ik 2 V [λ 2 + k 2 2 ]U -1 W 1 λ 2 + k 2 2 V -1 F 1 +V [λ 2 + k 2 2 ]U -1 F 2 . (10.62)
Defining the two matrices the relationship (10.59) becomes

C =     -ik 2 V [λ 2 + k 2 2 ]U -1 W 1 λ 2 + k 2 2 V -1 0 0 -ik 2 V [λ 2 + k 2 2 ]U -1 W 1 λ 2 + k 2 2 V -1     (10.63) and D = V [λ 2 + k 2 2 ]U -1 0 0 V [λ 2 + k 2 2 ]U -1 , ( 10
A F 1 (0) F 1 (h) = BC F 1 (0) F 1 (h) + BD F 2 (0) F 2 (h) . (10.65)
Finally, according to the definition (10.57), the R matrix is given by

R = 1 A -BC BD . (10.66)
It is stressed that the R matrix is numerically stable. Indeed, exponential functions always have arguments such that they decrease and, when inverted, they are always added to wellconditioned functions. For example, in the matrix [A -BC], one can check that the diagonal blocs are well-conditioned since they do not contain any exponential function, while the offdiagonal blocs decrease exponentially. When inverted, this matrix will have the same behaviour with well-conditioned diagonal blocs and exponentially decreasing off-diagonal blocs.

R matrix for a stack of lamellar layers

Here, a system made of two lamellar layers is considered. The first layer is located between the planes x 3 = -h 1 and x 3 = 0, and the second layer between the planes x 3 = 0 and x 3 = h 2 . Let R 1 and R 2 be the R matrices associated with these layers:

F 1 (-h 1 ) F 1 (0) = R 1 F 2 (-h 1 ) F 2 (0) , F 1 (0) F 1 (h 2 ) = R 2 F 2 (0) F 2 (h 2 ) . (10.67)
Then, the R matrix associated with the stack of the two layers is determined by eliminating the components F 1 (0) and F 2 (0) in the equations above. Denoting by R 1,i j and R 2,i j (i, j = 1, 2) the blocs of R 1 and R 2 ,

R 1 = R 1,11 R 1,12 R 1,21 R 1,22 , R 2 = R 2,11 R 2,12 R 2,21 R 2,22 , (10.68)
the expression of R is given by

R =     R 1,11 -R 1,12 1 R 1,22 -R 2,11 R 1,21 R 1,12 1 R 1,22 -R 2,11 R 2,12 -R 2,21 1 R 1,22 -R 2,11 R 1,21 R 2,22 -R 2,21 1 R 1,22 -R 2,11 R 2,12     .
(10.69) Again, one can check that the algorithm is stable since the only inverted blocs are the diagonal ones, which are well-conditioned.

Numerical application

A simple numerical example is considered to put the exact modal method to the test. The structure is made of a set of rectangular rods with dielectric constant ε 1,1 /ε 0 = 12.96 (corresponding to the index 3.6 of Si at optical wavelengths), width w 1,1 = 0.28d and height

h = d/(2 √ 2)
, where d is the spatial period of the grating (see figure 10.4). This lamellar grating is illuminated by a plane wave with an incident angle θ i = 45 • . The oscillating frequency is ω = 2π/(d √ ε 0 µ 0 ), which corresponds to a wavelength equal to the spatial period d. The efficiency diffracted in the order zero, i.e. at the reflected angle θ r = θ i = 45 • is calculated for both s and p-polarizations which correspond respectively to the electric and magnetic fields reduced to a single component along the invariance axis x 2 . Each component of the electromagnetic field is described by a finite number (2n + 1) of exact modes. Reflected efficiency in the order zero for different values of the number of exact modes (2n + 1) is represented on figure 10.5. These curves show that these efficiencies differ from their converged value with less than one percent from (2n + 1) = 7 in s-polarization and (2n + 1) = 5 in p-polarization (the converged values are respectively 0.7323 and 0.9487 in s and p polarizations). This convergence is found to be faster than in the case of the modal method with Fourier basis [START_REF] Li | New formulation of the Fourier modal method for crossed surface-relief gratings[END_REF] where an error smaller than one percent is obtained from (2n + 1) = 17 and (2n + 1) = 27 in s and p polarizations respectively (the method [START_REF] Li | New formulation of the Fourier modal method for crossed surface-relief gratings[END_REF] contains all the techniques to improve the convergence of the truncated Fourier series [START_REF] Li | Use of Fourier series in the analysis of discontinuous periodic structures[END_REF]). This improvement of the convergence resulting from the use of the exact modes becomes a significant advantage when three dimensional woodpile structures are considered [START_REF] Gralak | Theoretical study of photonic band gaps in woodpile crystals[END_REF] (the total number of modes (2n + 1) 2 can be reduced by a factor of 10).

x 3 x 1 ε 1,1 w 1,1 w 1,2 d ε 0 , µ 0 ε 0 , µ 0 θ i θ r

Appendix. Calculation of the exact modes and eigenvalues

It is shown here how to determine exactly the eigenvalues and the eigenfunctions of the operator K associated with a lamellar layer in a very general case. An analogous reasoning provides the ones of the operator associated with the others lamellar layers.

From the expression (10.50), every eigenvalue Λ n of the operator K is either an eigenvalue of

K ε 1 or K µ 1 . So, it is sufficient to determine the set of eigenvalues {Λ ν 1 ,n | n ∈ N } associated with the set of eigenfunctions {φ ν 1 ,n | n ∈ N } of the scalar operator K ν 1 , with ν 1 = ε 1 or ν 1 = µ 1 .

The equation satisfied by the exact eigenvalues

From the expression (10.50), the operator

K ν 1 is the sum of ω 2 ε 1 µ 1 + ∂ 1 ν -1 1 ∂ 1 ν 1 and ∂ 2 2
: the first part is an operator of the single variable x 1 and the second part is an operator of the single variable x 2 . Thus, we can perform a variable separation: every eigenfunction of L ν 1 can be written

φ ν 1 ,n (x 1 , x 2 ) = φ (1)
n 1 (x 1 ) φ (2)
n 2 (x 2 ) n 1 , n 2 ∈ N , (10.70) where φ (1)
n 1 and φ

n 2 are respectively eigenfunctions of the first and second operators which constitute L ν 1 .

In the case of a lamellar grating which is invariant in the direction x 2 , this direction of invariance is considered using the Fourier decomposition (10.5). Thus the eigenfunction of

∂ 2 2 is just φ (2) n 2 (x 2 ) = exp[ik 2 x 2 ] k 2 ∈ R , (10.71)
and the integer n 2 plays no role (integer n will be simply n 1 ). In the case of woodpile crystals, it is easy to verify that the plane-wave φ

n 2 (x 2 ) = exp{i[k 2 + 2π p(n 2 )/d 2 ]x 2 /d 2 } p(n 2 ) ∈ Z (10.72) (2)
is an eigenfunction of the operator ∂ 2 2 and satisfies the partial Bloch boundary condition (10.8) adapted for the variable

x 2 . Let Λ (2)
n 2 be the associated eigenvalue. Then, from (10.71,10.72), Λ

n 2 = -[k 2 + 2πq(n 2 )/d 2 ] 2 . (2)
(10.73)

Note that, for lamellar gratings and eigenfunctions (10.71), the integer q(n 2 ) is set to zero. The x 1 -dependency of the eigenfunction (10.70) is determined using the usual transfer matrix [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF][START_REF] Figotin | Localized electromagnetic waves in a layered periodic dielectric medium with a defect[END_REF][START_REF] Felbacq | Wave propagation in one-dimensional photonic crystals[END_REF]. Let λ (1) n 1 be the eigenvalue associated with φ (1)

n 1 : ω 2 ε 1 µ 1 + ∂ 1 ν -1 1 ∂ 1 ν 1 φ (1)
n 1 = Λ (1)
n 1 φ (1)
n 1 . (10.74)
In order to obtain a set of first order differential equations, the following column vector is introduced

F n 1 = ν 1 φ (1) n 1 ν -1 1 ∂ 1 ν 1 φ (1) n 1 .
(10.75) Note that, from equation (10.74), the two components of this vector are continuous functions. Now, suppose that the unit cell of the considered lamellar layer is made of J rods of width w 1, j , permittivity ε 1, j and permeability µ 1, j , j = 1, 2, • • • , J (figure 10.6): we denote by ν 1, j the value

x 3 x 1 ν 1,1 ν 1,2 ν 1,3 ν 1 = ε 1 , µ 1 w 1,1 w 1,2 w 1,3 d Figure 10
.6: A layer made of three rods per unit cell (J = 3): the three rods have width w 1, j , permittivity ε 1, j and permeability µ 1, j , j = 1, 2, 3.

of the function ν 1 in the rod j, j = 1, 2, • • • , J. Then, from equation (10.74), the vector (10.75) satisfies [START_REF] Felbacq | Wave propagation in one-dimensional photonic crystals[END_REF] F n 1 (d) = T 1 Λ

(1)

n 1 F n 1 (0) , (10.76)
where

T 1 (Λ) = T 1,J (Λ) T 1,J-1 (Λ) • • • T 1,1 (Λ) , (10.77)
T 1, j (Λ) = P 1, j (Λ, w 1, j ) , (10.78) P 1, j (Λ, w) = cos(β 1, j w) ν 1, j β -1 1, j sin(β 1, j w) -ν -1 1, j β 1, j sin(β 1, j w) cos(β 1, j w) , (10.79)

β 1, j = ω 2 ε 1, j µ 1, j -Λ j = 1, 2, • • • , J . (10.80)
Note that the four elements of each matrix T 1, j only depend on β 2 1, j : the expression (10.79) is independent of the definition of the square root (10.80). In addition to (10.76), the vector (10.75) has to satisfy the partial Bloch boundary condition (10.8) for the variable x 1 :

F n 1 (d) = exp[ik 1 d]F n 1 (0) . (10.81)
The combination of (10.76) and (10.81) implies that exp[ik 1 d] is an eigenvalue of the matrix T 1 Λ

(1)

n 1 : the equation det T 1 Λ (1)
n 1 . This last equation can be simplified using the fact that det T 1 = 1 (since, from (10.78), det T 1, j = 1, j = 1, 2, • • • , J): if exp[ik 1 d] is an eigenvalue of T 1 , then exp[-ik 1 d] is also. Thus, the equation (10.82) is equivalent to

trT 1 Λ (1) n 1 -2 cos[k 1 d] = 0 , (10.83)
where trT 1 is the trace of matrix T 1 . Once the eigenvalues Λ

n 1 are determined from (10.83), the associated eigenvectors φ (1) n 1 are also obtained using the transfer matrix [START_REF] Figotin | Localized electromagnetic waves in a layered periodic dielectric medium with a defect[END_REF]: firstly, the eigenvector F n 1 (0) in C 2 (associated with the eigenvalue exp[ik 1 d]) of the matrix T 1 Λ

(1) n 1 is determined; secondly, the expression of φ

(1)
n 1 in the rod j can be deduced from

F n 1 (x 1 ) = P 1, j Λ (1) n 1 , x 1 -x 1, j F n 1 (x 1, j-1 ) , (10.84)
where

x 1,0 = 0 , x 1, j = j ∑ q=1 w 1,q j = 1, 2, • • • , J . (10.85)
Finally the eigenvalues of the operator K ν 1 are

λ ν 1 ,n = λ (1)
n 1 + λ (2) n 2 , (10.86)
whose the two parts are respectively given by (10.83) and (10.73), and the expression of associated eigenvectors is (10.70) whose the two parts are respectively given by (10.84) and (10.72).

Concerning the functions of the operators U = σ 1 V and W = σ -1 1 ∂ 1 σ 1 V used in section 10.4 [see equation (10.56)], they are equal to the functions

ν 1 φ (1) n 1 (x 1 ) φ (2) n 2 (x 2 ) , ν -1 1 ∂ 1 ν 1 φ (1)
n 1 (x 1 ) φ (2)
n 2 (x 2 ) , (10.87)
where n 1 and n 2 are in N, the expression of ν 1 φ

(1)

n 1 and ν -1 1 ∂ 1 ν 1 φ (1)
n 1 in the rod j can be deduced from (10.75,10.84) and the expression of φ (2) n 2 is given by (10.72).

Real eigenvalues

Here, we suppose that the permittivity and permeability are real positive functions:

ε 1 (x 1 ) ∈ R , ε + > ε 1 (x 1 ) > 0 ; µ 1 (x 1 ) ∈ R , µ + > µ 1 (x 1 ) > 0 . (10.88)
Under these conditions, the operator K ν 1 is selfadjoint and its eigenvalues are real when the following inner product is used:

(φ , ψ) -→ 1 d d 0 ν 1 (x 1 )φ (x 1 )ψ(x 1 ) dx 1 . (10.89)
The only difficulty in the numerical determination of the eigenvalues (10.86) is to find the real numbers λ

n 1 which satisfy the transcendental equation (10.83). Since the numbers Λ (1)

n 1 are eigenvalues of the operator ω 2 ε 1 µ 1 + ∂ 1 ν -1 1 ∂ 1 ν 1 ≤ ω 2 ε + µ + ,
these numbers are on the semi-axis (-∞, ω 2 ε + µ + ]. This property makes their numerical determination easier and provides a way to number them:

ω 2 ε + µ + ≥ Λ ν 1 ,1 ≥ Λ ν 1 ,2 • • • ≥ Λ ν 1 ,n ≥ • • • (10.90)
However, two difficulties can occur in this numerical determination. We give herein the solutions we have adopted.

The first difficulty comes from the possibility for two consecutive numbers Λ The expression of the derivative of the matrix T 1 can be deduced from (10.77,10.78):

dT 1 dΛ = dT 1,J dΛ T 1,J-1 • • • T 1,1 + T 1,J dT 1,J-1 dΛ • • • T 1,1 + • • • + T 1,J T 1,J-1 • • • dT 1,1 dΛ , (10.92)
where, for j = 1, 2, • • • , J, the derivative of matrices

dT 1, j dΛ = 1 2 a 1, j b 1, j c 1, j d 1, j (10.93) is given by a 1, j = w 1, j β -1 1, j sin[β 1, j w 1, j ] , b 1, j = ν 1, j β -3 1, j sin[β 1, j w 1, j ] -ν 1, j w 1, j β -2 1, j cos[β 1, j w 1, j ] , c 1, j = ν -1 1, j β -1 1, j sin[β 1, j w 1, j ] + ν -1 1, j w 1, j cos[β 1, j w 1, j ] , d 1, j = w 1, j β -1 1, j sin[β 1, j w 1, j ] . (10.94)
The second difficulty comes from the possibility of numerical instabilities in the expressions (10.79,10.93) since the numbers β 1, j (10.80) can have non-vanishing imaginary part. A solution is to multiply the four coefficients of matrices T 1, j and their derivative (10.93) by the number

N j = exp -| Im(β 1, j )| w 1, j j = 1, 2, • • • , J , (10.95)
and the term 2 cos[k 1 d] wich appears in (10.83) by the product

N = N J N J-1 • • • N 1 . ( 10
.96)

Complex eigenvalues

Here, the permittivity and permeability can take any complex value: ν 1, j is in C, where ν 1 = ε 1 , µ 1 and j = 1, 2, • • • , J. The operator K ν 1 is not selfadjoint and then, its eigenvalues are, in general, in the complex plane. The determination of these complex eigenvalues λ

n 1 which satisfy the equation (10.83) has been intensively studied using different methods [START_REF] Botten | Highly conducting lamellar diffraction grating[END_REF][START_REF] Tayeb | On the numerical study of deep conducting lamellar diffraction grating[END_REF][START_REF] Sandström | Lossy multistep lamellar gratings in conical diffraction mountings: an exact eigenfunction solution[END_REF].

We present here a method similar to the one presented in [START_REF] Tayeb | On the numerical study of deep conducting lamellar diffraction grating[END_REF]: the complex eigenvalues are deduced from the real eigenvalues by an analytic continuation. However, our method differs from the one presented in [START_REF] Tayeb | On the numerical study of deep conducting lamellar diffraction grating[END_REF] since we make varying the phase of the numbers ν 1, j instead of their imaginary part. We think that it is better to make varying the phase since, from that we have observed, it leaves invariant the generalization to the complex case

Re(Λ ν 1 ,1 ) ≥ Re(Λ ν 1 ,2 ) • • • ≥ Re(Λ ν 1 ,p ) • • • (10.97)
of the numbering used when the eigenvalues are real (10.90).

We define for all t in [0, 1] the functions

ν 1, j (t) = |ν 1, j | exp[it arg(ν 1, j )] , (10.98)
where arg(ν 1, j ) is the phase of the complex number ν 1, j , ν 1 = ε 1 , µ 1 and j = 1, 2, • • • , J. Substituting the numbers ν 1, j (where ν 1 = ε 1 , µ 1 ) for ν 1, j (t) in equations (10.77,10.78), we obtain the matrix T 1 (Λ,t). For each value of t, we define the numbers Λ

n 1 (t) which satisfy tr T 1 Λ

(1)

n 1 (t),t -2 cos[k 1 d] = 0 . (10.99)
Then, the numbers Λ

n 1 (1) are the desired complex eigenvalues Λ

p 1 and the numbers Λ

n 1 (0) are real eigenvalues which can be determined using the method presented in the previous section 10. [START_REF] Tayeb | On the numerical study of deep conducting lamellar diffraction grating[END_REF] where ∂ T 1 /∂ Λ is given by substituting the numbers ν 1, j for ν 1, j (t) in equations (10.94,10.93) and ∂ T 1 /∂t is determined similarly. Finally the obtained estimates of numbers Λ

n 1 (1) are used to initiate any of the classical methods for the numerical solution of equations [START_REF] Tayeb | On the numerical study of deep conducting lamellar diffraction grating[END_REF]. Then, one obtains the desired complex eigenvalues.

In order to eliminate the numerical instabilities, one has to multiply each matrix T 1, j and their derivatives by the numbers N j (10.95) as in the previous section 10.7.2.

Eigenfunctions

From (10.84), the expression of each eigenfunction φ (1) n 1 is given by the coefficients of the column vectors F n 1 (x 1, j ), j = 0, 1, • • • , J. On the numerical side, the only difficulty comes from the fact that numerical instabilities in the expression of the transfer matrices (10.78,10.79). A solution based on the R-matrix algorithm (or S-matrix) should consist in using the algorithm presented in [START_REF] Gralak | From scattering or impedance matrices to Bloch modes of photonic crystals[END_REF] to obtain the vector F n 1 (x 1,0 ) (and the vector F n 1 (x 1,J ) = exp[k 1 d]F n 1 (x 1,0 )) and then, the algorithm presented in [21, section V] to obtain the vectors F n 1 (x 1, j ), j = 1, 2, • • • , J -1. However, we propose to use another solution which benefits of the fact that we deal with 2 × 2 matrices.

We define the following complex coefficients: = T 1, j Λ

n 1 T 1,J-1 Λ (1) n 1 • • • T 1, j Λ (1)
(1) where, as in (10.104), the multiplication by the number N consists in multiplying each coefficient of matrices T 1, j Λ

n 1 T 1, j-1 Λ (1)
n 1 • • • T 1,1 Λ (1)
n 1 by the number N j . Finally these functions have to be normalized. From the definition (10.89) of the inner product, one has to compute (10.107)

This expression allows to eliminate the numerical instabilities which can occur from the exponential functions. Note that all the coefficients of matrices defined in sections 10.4 and 10.5 (matrices U, V and W ) can be also computed analytically in order to eliminate the numerical instabilities.

Figure 10 . 2 :

 102 Figure 10.2: A lamellar grating made of a single lamellar layer on a substrate.

Figure 10 . 3 :

 103 Figure 10.3: Integration in the complex plane of k 3 .

Figure 10 . 4 :

 104 Figure 10.4: The considered structure for the numerical example: a single layer made of rectangular rods.

Figure 10 . 5 :

 105 Figure 10.5: Efficiency in the zero order for s-polarization (left panel) and p-polarization (right panel).

n 1 -

 1 exp[ik 1 d] = 0 (10.82) determines the eigenvalues Λ

  be very close to each other. Our solution is to use an algorithm which determines the zeros of the function trT 1 (Λ) -2 cos[k 1 d] on the left side of equation (10.83) by taking into account this function together with its derivative with respect to Λ. If two numbers Λ (1) n 1 are very close one to each other, then the derivative is close to zero. Thus such algorithm needs to determine the function d dΛ trT 1 (Λ) -2 cos[k 1 d] = tr dT 1 dΛ (λ ) . (10.91)

1 F j 2 = 11 N

 1211 F n 1 (x 1, j ) j = 0, 1, • • • , J . (10.103) Since F n 1 (x 1,0) is an eigenvector of the matrix T 1 Λ(1)n 1 associated with the eigenvalue exp[k 1 d]), and T J 12 N are obtained by multiplying each coefficient of matrices T 1, j Λ(1) n 1 by the number N j . The coefficients F J 1 and F J 2 are deduced from (10.81,10.104) and then, one can obtain the other coefficients for j = 1, 2, • • • , J -1:

  φ

n 1 (x 1 ) 2 ν 1

 121 (x 1 ) dx 1(10.106) when the functions ε 1 and µ 1 have the property (10.88). In the general case (where ε and µ are complex valued functions), one has to use the formalism presented in [8, section 2.3]. It is possible to compute analytically the expression (10.106):

  Table of Contents: 10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.2 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.3 Continuation of the electromagnetic field . . . . . . . . . . . . . . . . . . . . . 10.3.1 Direct formulation: the transfer matrix . . . . . . . . . . . . . . . . . . 10.3.2 Rigorous derivation of the continuation procedure . . . . . . . . . . . . 10.4 Exact eigenmodes and eigenvalues method . . . . . . . . . . . . . . . . . . . . 10.5 Numerical algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.5.1 R matrix for a single lamellar layer . . . . . . . . . . . . . . . . . . . . 10.5.2 R matrix for a stack of lamellar layers . . . . . . . . . . . . . . . . . . 10.6 Numerical application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.7 Appendix. Calculation of the exact modes and eigenvalues . . . . . . . . . . .

10.7.1 The equation satisfied by the exact eigenvalues . . . . . . . . . . . . . 10.7.2 Real eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.7.3 Complex eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.7.4 Eigenfunctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  .64)

	10.12	Gratings: Theory and Numeric Applications, 2012

  .2. Assuming that Λ

	(1) n 1 (t) are continuous and differentiable functions of t, the complex
	numbers Λ (1) n 1 (1) can be estimated from the numbers Λ (1) n 1 (0) by a numerical integration [7] of
	d Λ dt (1) n 1	(t) = -	tr(∂ T 1 /∂ Λ) Λ (1) n 1 (t),t tr(∂ T 1 /∂t) Λ n 1 (t),t (1)	,	(10.100)


Stephanie Deboeuf

Eytan Katzav

Arezki Boudaoud

Daniel Bonn

Mohktar Adda-Bedia

A comparative study of crumpling and folding of thin sheets

Keywords: numbers: 46.25.-y, 46.32.+x, 46.70.-p, 62.20.-x

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

It is easy to verify that the maximum number of times one can fold a sheet of paper is only 6 or 7, which surprisingly is independent of the initial size of the sheet. Quantitatively, elasticity theory allows to write the relation between the compaction force and the number of times one can repeatedly fold a piece of paper in two. This follows from the scaling of the bending rigidity B with the thickness h of the folded sheet [START_REF] Landau | Theory of Elasticity[END_REF], B = Eh 3 /12(1 -ν 2 ), where E is the Young modulus and ν the Poisson ratio. For a sheet of initial size D 1 × D 2 folded along the direction D 2 , the compression energy, E c , injected in the system should be compared to the typical energy dissipated in the fold. One writes E c = F D 2 , where F is a characteristic compression force applied along the direction D 1 . Since most of the folded sheet remains flat and the region which is irreversibly deformed is straight along D 2 (i.e. its gaussian curvature equals to zero ; See Fig. 1), the energy dissipated in the fold is estimated from the elastic bending energy, E el , concentrated in a region of length D 2 and width h with a curvature 1/h [START_REF] Landau | Theory of Elasticity[END_REF]; this leads to E el ∼ BhD 2 /h 2 = BD 2 /h [START_REF] Dias | We can retrieve the energy of the fold by using the estimate of the energy of the crease given[END_REF]. The balance of these two energies leads to F 0 ∼ B/h for the elementary force needed to create a unique fold. When the sheet is folded n times repeatedly leading to the hierarchical creation of folds, its effective thickness and bending rigidity become h n → 2 n h and B n → 2 3n B assuming no slip between layers (this hypothesis becomes increasingly consistent for large n). Consequently, for a sheet folded n times, the energy balance gives

F (n) ∼ B n /h n ∼ F 0 2 2n . ( 1
)
Thus the force is independent of the initial size of the sheet and grows exponentially with the number of folding events n. The exponential dependence is the reason why one cannot fold a sheet indefinitely by hand or by applying a finite force; the elementary force F 0 is estimated with typical values E = 10 9 Pa, h = 10 -4 m leading to F 0 ≈ 1 N. Then, for n = 6, F becomes of the order of kiloNewtons, which is larger than the maximal force any person can exert and then sets a limit on the achievable number of successive foldings. Repeated folding in two is not the only possible way to fold; here we consider three basic regular processes (Fig. 1), which are prototypical foldings in various dimensionalities. In cases (a) and (b), the sheet is thought of as a 1d-like sheet, since it is folded along one direction only. However compaction of case (a) is not isotropic, contrary to case (b). The latter is seen as isotropic compaction within a disk, whereas the former is seen as unidirectional compaction inside an elongated rectangle such as in [3]. Finally, for case (c) the sheet is truly twodimensional and compacted in a sphere. The number of folded layers N = h n /h after n folding events and the related compaction ratio φ (defined as φ ≡ D/∆, where D and ∆ are the initial and final size of the sheet) then depends on the precise geometry and dimensionality of the compaction process. Folding a d-dimensional sheet in (d + 1)-dimensional space, N = 2 n , whereas for the 1d-1d case (a), N = n. Also, for foldings (a) and (b), N = φ, whereas in the 3d-like folding (c), N = φ 2 . Using similar arguments leading to Eq. ( 1), one finds a generic power law relation between the force F , the compaction ratio φ and the number of folded layers N : 1d compaction, case (a) :

F (N ) = F 0 N ∼ F 0 φ (2) 2d compaction, case (b) : F (N ) = F 0 N 2 ∼ F 0 φ 2 (3) 3d compaction, case (c) : F (N ) = F 0 N 2 ∼ F 0 φ 4 (4)
It is important to note that we described folds as the result of an irreversible process occurring in a small region of size h and characterized by a zero gaussian curvature. Consequently, energy scalings are different from those obtained for singular ridges [4,5] that are reversible. Moreover, our estimate for the energy dissipated in the fold should be taken as a lower bound because the contribution of vertices [6] and other possible length scales associated with plastic events are neglected. Our approach is inspired by crack modelling in the framework of linear elastic fracture mechanics where the process zone near the crack tip is neglected and the dissipation is estimated from its balance with the far-field elastic energy.

One may wonder to which extent our models for regular folding describe our crumpling measurements and whether crumpling process shall be viewed as arising from successive folding events. To compare the regular folding with the crumpling of paper, we first show experimentally that here the force also increases as a power law with the degree of compaction, with exponents in accordance with predictions dependant on the dimensionalities of the compaction process. For this purpose, a sheet of paper of characteristic size D, is placed into a rigid cylindrical cell, in which a piston connected to a force transducer compacts the sheet at constant speed (Fig. 2a). The materials used here (Kraft and regular printing papers) have been chosen for their low ductility. The experimental force-distance curves show a very strong increase of the force upon compaction (Fig. 2b), for both 'virgin' sheets (crumpled for the first time) and 'trained' sheets (crumpled for the second or third time). For the latter the force-compaction ratio curves turn out to be independent of the initial preparation of the sheet inside the cylinder within the experimental accuracy. Then, the measured curves for different types of paper of similar properties (thickness h = 10µm and Young modulus E 10 9 P a), different sheet sizes, cells and compaction speeds are all described by a power law: where α is a characteristic force scale and ∆ is the gap between the piston and the bottom of the cell (Fig. 2a). While the data range for ∆ is small, this behaviour is robust over all the 150 realizations. The statistical χ 2 test for goodness-of-fit confirms the relevance of the power law in comparison with other fits. The exponent β of the power law divergence is β ≈ 1.3 (Fig. 2c), a value between 1 and 2, those expected for ordered folding in 1d and 2d. We argue below that this is due to the anisotropy of the compaction process in our experiment. Effectively, compaction here is quasi 1d, since loading is applied mainly in one direction. However, the setup also allows for compaction in the perpendicular direction, which would rather be a 2d process. Moreover, the characteristic force scale α is independent of size D: α ≈ 2 N (Fig. 2d) which is of the same order of magnitude as the characteristic force, F 0 ≈ 1 N, calculated for the folding.

F (∆) = α ∆ D -β = αφ β , (5)
A second step towards understanding the analogy between crumpling and folding is to establish the relation between the degree of compaction and the number of folds for the experimental crumpled configurations. To achieve this, we characterize the geometry of the crum- pled paper, using an original approach that makes use of the properties of folds and facets in cross-sections of crumpled samples. Sheets of different paper types and sizes D are crumpled into hand-made balls at different degrees of compaction. A cross-section is obtained by cutting the crumpled ball in two with a slowly moving hot wire [7]. The overall size of the crumpled configuration, ∆, is defined as the largest diameter of the resulting cross-sectional area. In this cross-section, the number of paper layers is measured in two orthogonal directions passing through the center (Fig. 3a), and averaged to obtain the mean number of folded layers N in the crumpled configuration. By this method, we ensure that N is defined as in the folding model introduced above (N = h n /h). The number of folded layers N is described by a power two dependance on the degree of compaction N = (D/∆) 2 = φ 2 , with a prefactor equal to 1 (Fig. 3c). This important result is exactly the same as that observed for 3d folding (Fig. 1c) showing that the geometry of folding and crumpling is the same. A possible difference between the two situations is that while the repeated folding is a hierarchical process, this is not clear for crumpling. To investigate whether the crumpling is also hierarchical, we charac-terize the lengths of folds and facets in cross-sections of crumpled samples. For this purpose, the cut crumpled sheet is reopened carefully and thel uncrumpled pieces, with possibly several holes, are scanned (Fig. 3b). The edges of their boundaries and holes are detected automatically and broken down into segments delimited by kinks [8], by using a 'Split and Merge' algorithm [9] for the segmentation. The planar two-dimensional crosssection of the crumpled sheet bears information on the full three-dimensional crumpled configuration: the ensemble of segments samples the facets delimited by folds, so that one segment length is related to the characteristic size of the facet or equivalently to the characteristic distance between folds. We use this to asses the nature of the crumpling. To do so, we compare the distribution of lengths ρ( ) with a log-normal distribution and a distribution; the former characterizes a hierarchical process [4,8,10], whereas the latter accounts for random processes [4,8,10]. More precisely, a lognormal distribution describes a fragmentation process in which all pieces are broken successively into two parts, such that any new fragment is further broken into two pieces where the breaking point is uniformly distributed along the fragment [8,11,12]. In contrast, the Gamma distribution emerges from a fragmentation process where all the breaking points are uniformly distributed along the unbroken line, prior to the breaking that happens simultaneously for all points. Fig. 3d shows that both distributions reasonably well describe the rapid decay of the tail of the distribution, but the Gamma distribution seriously overestimates the probability density at small lengths. This originates from the fact that a hierarchical fragmentation process tends to generate less small fragments than a random one. A more rigorous test is done through the statistical χ 2 test for goodness-of-fit, which confirms that the log-normal describes better the data. We checked that this description is robust with respect to the chosen value of the threshold used in the segmentation procedure. The log-normal distribution accurately describes all the experimental data sets, so the crumpling is hierarchical rather than random. Earlier simulations [13] of crumpled sheets and experiments on unfolded sheets [14] found a similar agreement with a log-normal distribution.

The conclusion is that folding and crumpling are very similar in nature and the crumpling process shall be viewed as arising from successive folding events. For ordered folding, simple models allow for predictions of the relations between force, compaction ratio and number of folds. Surprisingly, these are found to capture the main properties of crumpling also, in particular the hierarchical structure of the folds and the power law relation between the force and the compaction ratio. The analogy with folding then allows to define the 'crumpliness exponent' β for various forms of crumpling process. Previous experiments and simulations in the literature have TABLE I. Results from the literature for the power-law variation of the force F ∝ (D/∆) β for compactions of an xdimensional object in (x + 1)-dimensional space. The crumpliness exponent β is measured and β is our theoretical prediction from the dimension x of the geometry of compaction. Topological constraints and the material properties may influence the value of the exponent. The first set corresponds to a mixture of cases (a) and (b) of Fig. 1 reported such exponents β for the power-law dependence of the force on the compaction ratio. They can be explained using our arguments, i.e., solely by considering the dimensionality of the compaction process, the topological constraints and the mechanical properties of the material (e.g. ductility). Table I summarizes various exponents β from the literature, detailed below.

For the first set of data, Matan et al. [17] used a compaction set-up similar to the one used here, and found an exponent of β ≈ 1.89. The aspect ratio of their cylinder (height/diameter) is much smaller than ours; we thus anticipate that the compaction is more 2d in nature, and hence expect a crumpliness exponent closer to 2 than in our experiment. The value of this exponent can again be understood as a compaction process lying between 1d and 2d. As our arguments are based on dimensionality, they allow to predict only bounds for this type of experiments. The simulations of [18] on compacted tethered membranes found a value of β ≈ 1.85. Except that loading is now biaxial, the compaction process is in fact similiar to that of case (a) since the "height" fluctuations of the membrane are small. If the two directions were independent we recover N ∝ φ 2 and thus force F = N F 0 ∝ φ 2 . However, folding in one direction is inhibited by folding in the other direction. This effect will decrease the total number of folds leading to a crumpliness exponent β 2. For the second set, the analogy with case (b) is complete and both experimental [20], theoretical [19] and numerical [20,21] results in the elastic regime are in agreement with the prediction of a crumpliness exponent β = 2.

The third set deals with experiments and simulations of 2d sheets crumpled inside 3d spheres. The linearly elastic sheet [13] is a perfect example of case (c) for which the crumpliness exponent β = 4, in agreement with the simulations. The aluminum foil is ductile with plastic deformations [22] and the phantom sheet can cross itself [13]. These are more complicated cases; however we can estimate the exponent β. For the aluminum foil, one modifies the estimate of the elastic energy E el of a folded sheet because of ductility. Most of the folded sheet remains flat but the elastic energy is now concentrated in a region of length D and width 1/κ c with a curvature κ c , which is a material constant: the curvature scale at which the material yields. The balance of the compaction energy and the elasto-plastic energy then leads to F ∼ B n κ c ∼ Eh 3 n κ c . As N = φ, one finds F (N ) ∼ N 3 F 0 ∼ F 0 φ 6 leading to the value β = 6, which is in fair agreement with experimental results [22]. For the phantom sheets, the absence of steric interactions implies that F (N ) ∼ N F 0 which is similar to the 1d case. The number of folded layers N is then related to the compaction degree through N ∼ V /V f , where V f is the average volume occupied by the sheet. For high compactions, it is known that

V f ∼ R d f
g where R g is the radius of gyration and d f is the fractal dimension [23]. Thus we find F ∝ N ∼ φ β , with β = d f 2.5 [23,24].

Finally, these arguments allow to explain why a wastebucket fills up so quickly when waste paper is crumpled into a ball. Using the equivalence between crumpling and folding, the wasted volume ∆V /V can be estimated from the folded case. In the 3d case, this is given by

∆V V ∆ 3 -N h∆ 2 ∆ 3 = 1 - h D N 3/2 . (6)
For N = 2 6 , h = 10 -4 m and D = 0.2m, one has ∆V /V ≈ 75%, which is an excellent estimate for the experimental observation that crumpling is a very inefficient compaction process. In conclusion, the observations presented here demonstrate a non-trivial relation between the force of compaction and the geometry of the crumpled configuration. A potential application of this result would be to invert this problem, and deduce the force through analyzing cross sections of crumpled sheets. Since our arguments are generic, they should hold at the nano-scale and could provide a simple framework to understand crumpled graphene structures, such as graphene-based supercapacitors [15,16].

FIG. 1 .

 1 FIG. 1. Hierarchical folding of a sheet in different dimensionalities. (a) a 1d-like sheet folded in 1d; (b) a 1d-like sheet folded in 2d; (c) a 2d-like sheet folded in 3d. The three types of folding processes are referred to as 1d, 2d and 3d compactions.

  FIG. 2. (a) Setup used to measure the force F (∆) during crumpling. (b) Typical force-distance F (∆) and forcecompaction ratio F (∆/D) curves on a log-log scale and shifted with respect to each other for clarity in inset, with Kraft paper of size D = 15cm in a cell of diameter 6cm, for the 1 st (blue circles), 2 nd (green squares) and 3 rd (red diamonds) crumpling rounds. The lines are the fits to a power law (Eq. 5): (α, β) are (3.74 N, 1.1), (1.90 N, 1.4) and (1.43 N, 1.5) for the 1 st , 2 nd and 3 rd crumpling rounds respectively. (c) Exponent β of the power law fit as a function of the paper size D and its probability distribution function in inset for the 2 nd (green) and 3 rd (red) crumpling rounds. (d) Characteristic force scale α of the power law crumpling force as a function of D and its probability distribution function in inset. For (c) and (d) many experimental realizations were averaged for constant values of D.

FIG. 3 .

 3 FIG. 3. (a) Picture of a crumpled cross-section and two orthogonal directions used to extract the number of folded layers N . (b) Picture of a piece of an uncrumpled cross-section and its segmented edge. Scale bars are 10mm. (c) Number of folded layers N as a function of compaction ratio D/∆ for crumpled balls of Kraft paper (circles) and regular printing paper (triangles). The line is the curve N = (D/∆) 2 . (d) Probability density function ρ( ) of the length of segments on a log-log scale for a Kraft paper sheet crumpled in a ball of diameter ∆ ≈ 80mm and compaction ratio φ ≈ 6. The continuous and dashed lines are the log-normal and Gamma distributions respectively, with the same mean and variance as the data. The error bars δ and δρ are the bin width δ and the standard deviation δρ = ρ/ √ n of the histograms n( ).


Samir F Matar
email: matar@icmcb-bordeaux.cnrs.fr

First principles study of the electronic and magnetic structures and bonding properties of UCoC 2 ternary, characteristic of C-C units

Keywords: Uranium cobalt di-carbide, Electronic Structure, EOS, Bonding

The electronic structure of UCoC 2 , a di-carbide with the C-C units is examined from ab initio with an assessment of the properties of chemical bonding. The energy-volume equation of state shows large anisotropy effects due to C-C alignment along tetragonal c-axis leading to high linear incompressibility. Relevant features of selective bonding of uranium and cobalt with carbon at two different Wyckoff sites and strong C-C interactions are remarkable. The vibrational frequencies for C…C stretching modes indicate closer behavior to aliphatic C-C rather than C=C double bond. A ferromagnetic ground state is proposed from the calculations.

Introduction

Uranium forms a series of binary compounds with carbon and silicon. The equiatomic UC with rock-salt structure contains isolated C atoms and no C-C bonds, while U 2 C 3 and UC 2 have C 2 pairs as carbon substructure. Similarly the higher congener silicon U 3 Si has isolated silicon atoms while U 3 Si 2 has Si 2 pairs [1]. A rare example of mixed uranium silicide carbide is U 3 Si 2 C 2 [2,3] which contains SiC units. Based on computations within the density functional theory (DFT) [4] the chemical bonding was identified to be specific to uranium sites with Si and C with major Si-C bonding [5]. At another level of theory Li and Hoffmann [6] discussed the chemical role of C…C units in several binary and ternary metal carbides with extended Hückel methodology, including title compound. X-ray diffraction studies of UCoC 2 by Gerss, W. Jeitschko led to assign a new structural type derived from PbFCl structure [7]. In continuation of our ab initio investigations of the electronic and magnetic structures and chemical bonding of uranium based compounds within DFT [5,[8][9][10][11], we address herein these properties for UCoC 2 with stressing the effect of the structural anisotropy brought by C…C pairs on site selective chemical bonding with metal atoms U and Co and on the mechanical properties. Further from energy criteria a ferromagnetic ground state is proposed.

Computational methodologies

Two computational methods within the DFT were used in a complementary manner.

The Vienna ab initio simulation package (VASP) code [12,13] allows geometry optimization and total energy calculations. For this we use the projector augmented wave (PAW) method [13,14], built within the generalized gradient approximation (GGA) scheme following Perdew, Burke and Ernzerhof (PBE) [15]. Also preliminary calculations with local density approximation LDA [16] led to an underestimated volume versus the experiment, i.e. 84 Å 3 versus 89.4 Å 3 [7]. The conjugate-gradient algorithm [START_REF] Press | Numerical Recipes[END_REF] is used in this computational scheme to relax the atoms. The tetrahedron method with Blöchl corrections [14] as well as a Methfessel-Paxton [START_REF] Methfessel | [END_REF] scheme were applied for both geometry relaxation and total energy calculations. Brillouin-zone (BZ) integrals were approximated using the special k-point sampling. The optimization of the structural parameters was performed until the forces on the atoms were less than 0.02 eV/Å and all stress components less than 0.003 eV/Å 3 . The calculations are converged at an energy cut-off of 400 eV for the plane-wave basis set with respect to the k-point integration with a starting mesh of 444 up to 888 for best convergence and relaxation to zero strains. The calculations are scalar relativistic.

In the context of the electronic structure investigations of actinides within the DFT, it may be argued that the correlations accounted for within its basic functionals (LDA and GGA) may not be sufficient enough to account for stronger correlations characterizing the 5f states. Then LDA+U or GGA+U method can be used to reproduce experimental findings. U stands for the Coulomb on-site repulsive parameter added to enhance correlation. Beside plain GGA we used such test calculations following Dudarev et al. [19], with U =4 eV value which was shown to better approach experiment in U 3 Si 2 C 2 [6].

Then all-electron calculations with the GGA were carried out for a full description of the electronic structure and the properties of chemical bonding, using full potential scalar-relativistic augmented spherical wave (ASW) method [20,[START_REF] Eyert | The Augmented Spherical Wave Method -A Comprehensive Treatment[END_REF]. In the minimal ASW basis set, we chose the outermost shells to represent the valence states and the matrix elements were constructed using partial waves up to l max +1=4 for U, l max +1 = 3 for Co and l max +1 = 2 for C. Self-consistency was achieved when charge transfers and energy changes between two successive cycles were below 10 -8 and 10 -6 eV, respectively. BZ integrations were performed using the linear tetrahedron method within the irreducible wedge. In order to optimize the basis set, additional augmented spherical waves are placed at carefully selected interstitial sites (IS). Besides the site projected density of states, we discuss qualitatively the pair interactions based on the overlap population analysis with the crystal orbital overlap population (COOP) [22]. In the plots, positive, negative, and zero COOP indicate bonding, anti-bonding, and nonbonding interactions, respectively.

Crystal chemistry

The crystal structure of UCoC 2 , space group P4/nmm, is sketched in Fig. 1 and described in Table 1. The unit cell can be considered as a stacking of Co, U and C…C along the tetragonal c-axis. It can be considered as an anti-PbFCl-type because Co occupies the F position at the origin (2a) whereas U is at Pb (2c) with 0 ½ z.

Carbon at chlorine sites (2c) is found within two sites, C1 and C2. It can be noted that this structure is different from actinide based homologue ThCoC 2 which crystallizes in the base centered monoclinic CeNiC 2 -type, indentified for a large number of compounds [23].

The C…C pairs (C1-C2, Fig. 1) have short distance of ~1.47 Å which is between C-C single bond of ~1.5 Å and C=C double bond (1.34 Å). Also the other short distances shown in Table 1 are those of U with C2 and Co with C1. Then significant bonding between them on one hand and between C1 and C2 on the other hand can be expected as it is discussed in section 5.2.

Results and discussions

Geometry optimization and equation of state (EOS)

Table 1 provides the experimental and calculated crystal structure parameters. With GGA the volume is found slightly below experimental value, in better agreement than with preliminary PAW-LDA calculations which severely underestimated it. The z parameters of U, C1 and C2 atoms at (2c) position are in good agreement with experiment as well as the interatomic distances which are found only slightly smaller.

The 3 rd column showing the GGA+U calculations lead to better agreement with experiment for the volume and for z(U) but to significant deviations of the z(C1) and z(C2) parameters as well. The energy is highly raised due to the U repulsive character.

Then it may be suggested that the changes are not significant of improvements on all structure properties as it was the case of U 3 Si 2 C 2 [5] and the following results are from calculations using plain PAW-GGA method.

For a better assessment of the results we derive the equilibrium zero pressure parameters from the energy-volume (E,V) equation of state (EOS) with (E,V) set of calculations around minima found from geometry optimization. The underlying physics of this procedure is that the calculated total energy corresponds to the cohesion within the crystal in as far as the solution of the Kohn-Sham DFT equations gives the energy with respect to infinitely separated electrons and nuclei. But the zero of energy depends on the choice of the potentials, then energy becomes arbitrary through its shifting, not scaling. However the energy derivatives as well as the EOS remain unaltered. For this reason one needs to establish EOS from which the fit parameters are extracted for an assessment of the equilibrium values. The resulting values are plotted in Fig. 2. They have a quadratic variation which can be fitted with an energy-volume Birch EOS up to the 3 rd order [24]: 3 , where E o , V o , B o and B ' are the equilibrium energy, the volume, the bulk modulus and its pressure derivative, respectively. The fit curves reproduce the trends of the geometry optimization for the energies and the volumes and provide the results given in the insert.

E(V) = E o (V o )+[9/8]V o B o [([(V o )/V]) [2/3] -1] 2 +[9/16]B o (B ' -4)V o [([(V o )/V]) [ 2/3] -1]
It can be noted that the volume is now closer to experiment. The corresponding zero pressure bulk modulus of B 0 = 260 GPa positions the hardness of the compound within the range of oxides. However in view of the anisotropic structure (Fig. 1) with C-C alignments we further investigated the linear compressibility along the c-axis. The obtained E,V curve shows a steeper change of the energy with volume leading to different fit parameters with the Birch EOS given above: the bulk modulus is strongly enhanced, B 0 = 560 GPa, with respect to isotropic volume compression while the equilibrium volume V 0 and energy E 0 show little change. This result may be assigned to the incompressibility of aligned C-C substructures with U. B' values are within range of usually encountered magnitudes [25].

Elastic stiffness constants

In tetragonal symmetry there are six independent elastic stiffness constants C 11 , C 33 , C 44 , C 66 , C 12 and C 13 . Most encountered compounds are polycrystalline where monocrystalline grains are randomly oriented so that on a large scale, such materials can be considered as statistically isotropic. They are then completely described by the bulk modulus B and the shear modulus G, which may be obtained by averaging over the single-crystal elastic constants. The most widely used averaging method is Voigt's based on a uniform strain (cf. [START_REF] Matar | Matériaux ultra-durs: Concepts et Modélisations[END_REF] for a review with therein refs.). The calculated elastic constants in UCoC 2 are (in units of GPa): The shear modulus which defines the rigidity of the material is significantly lower than the bulk modulus and the ratio G/B which is an indicator of brittleness/ductility has a value of 0.7. This defines a brittle material as for instance Ir metal contrary to neighboring transition metals Ag, Pt, au with G/B ratios in the range {0.4 -0.2} [27].

Analysis of the charge trends

Analyzing the charge density issued from the self consistent calculations can be done using the AIM (atoms in molecules theory) approach [28] developed by Bader who devised an intuitive way of splitting molecules into atoms as based purely on the electronic charge density. Typically in chemical systems, the charge density reaches a minimum between atoms and this is a natural region to separate them from each other.

Such an analysis does not constitute a tool for evaluating absolute ionizations. The We further assess the nature of the C…C bond with the calculation of the vibrational frequencies of the whole crystal system calculated through the Hessian matrix using finite differences with VASP package with parameter IBRION = 6 in INCAR control file (cf. web ref. in [13]). The results are such that at low frequencies (~150 cm -1 ) the lattice modes are found involving the displacements of the whole structure in a rigid manner and at the highest frequencies the C…C antisymmetric and symmetric stretching modes are obtained with (C…C) = 958 cm -1 and symmetric (C…C) = 1000 cm -1 respectively. Screen snapshots from the relevant stretching animations are shown in Fig. 3 with arrows indicating the directions. The highest value comes closer to infra red (IR) frequency for aliphatic C-C with (IR) = 1200 cm -1 than ethylene like C=C with (IR) = 1650 cm -1 . Note that a better agreement cannot be expected in as far as the vibration frequencies were calculated in the solid state whereas the text books refer to C-C and C=C in molecules.

Electronic structure and chemical bonding

Electronic density of states.

Calculations used the experimental parameters with the full potential scalar relativistic ASW method in which chemical bonding qualitative descriptions implementations were done by V. Eyert following different schemes [START_REF] Eyert | The Augmented Spherical Wave Method -A Comprehensive Treatment[END_REF].

Firstly spin degenerate (non spin polarized NSP) calculations with total spins are done. At self consistent convergence the charge transfer follows the trends observed above with additional charge residues (corresponding to less than 0.15 electrons) from the atomic spheres to IS. The site projected densities of states (PDOS) are shown in Fig. 4a. In both panels, the zero energy along the x-axis is considered with respect to the

Fermi level E F in as far as the ternary compound is metallic due to the PDOS contributions from all constituents, especially the lower part of the uranium 5f states.

Due to their low filling U(5f) are centered above E F . The opposite is observed for largely filled d states of Co which are centered below E F . Low energy lying carbon s at ~15 eV are followed by p PDOS around -10, -5 and from -4 up to E F with similar shapes to U and Co PDOS, thus ensuring for the chemical bonding with them.

While cobalt shows a small magnitude PDOS at E F , the lower part of U(5f) crosses E F at a large PDOS. This signals instability of the compound for uranium f states in such spin degenerate magnetic configuration within the Stoner mean field theory of band ferromagnetism [START_REF] Mohn | Magnetism in the solid state -An introduction[END_REF] which can be applied to address the tendency for spin polarization (not the long range magnetic order). The total energy of the spin system results from the exchange and kinetic energies. Referring the total energy to the non-magnetic state (NSP), this is expressed as: E = constant {1-I n(E F )}. In this expression, I (eV) is the Stoner integral and n(E F ) (1/eV) is the PDOS value for a given state -here f-at the Fermi level in the non-magnetic state (Fig. 4a). If the unit-less Stoner product I.n(E F ) is larger than 1, energy is lowered and the system stabilizes in a magnetically ordered configuration. Then the product I.n(E F ) provides a criterion for the stability of the spin system. From quantum theoretical calculations [30] the value of I{U(5f)} = 0.51 eV was derived. With n(E F ) value of 3.8 eV -1 the calculated I.n(E F ) is 1.94 and the Stoner criterion 1-I n(E F ) of -0.94 is negative, leading to energy lowering upon the onset of magnetization. Spin polarized (SP) calculations actually lead to a lowering of the total energy by 2.1 eV and to a total magnetization 1.3  B /formula unit in a ferromagnetic configuration. The major contribution to total magnetization arises from uranium with M(U)= 1.24  B while M(Co) = 0.06  B . Such magnetic polarization was suggested [5,7] without further specification. The effect on the PDOS is shown at Fig. 4b. In spin polarized calculations with DFT-GGA the magnetic moment arises from the rigid band shift between majority  spin populations and minority  spins ones through magnetic exchange. This is observed mainly for uranium f-PDOS around E F while the other constituents PDOS exhibit small (Co) or no (C1, C2) energy shift.

However the spin only value obtained from scalar relativistic calculations should be corrected by the spin orbit coupling SOC with the orbital moment which amounts to 2.7

in uranium [11]. This obeys Hund's 3 rd rule whereby the actual quantum number is J with J = |L-S| for less than half filled f subshell and J = |L+S| for more than half filled subshell. The first condition corresponds to uranium and the resulting moment is 2.4  B .

This should be confronted with experiment when available.

For the search of the actual magnetic ground state, an antiferromagnetic configuration was tested with dispatching the crystal into two magnetic subcells, one of them for  SPIN and the other for  SPIN. The result is zero magnetization as expected, and a raise of the energy by +0.5 eV/UCoC 2 unit. Thus the magnetic ground state is predicted to be ferromagnetic. Lastly we note that the magnitude of magnetization only slightly changed with GGA+U calculations (U = 4 eV) as in U 3 Si 2 C 2 [6]. This suggests band magnetism behavior of UCoC 2 ternary, due to the presence a transition metal of the first period T, (here Co) besides uranium like in UT 2 X 2 (X = p element) [9,10].

Chemical bonding.

Fig. 5 shows the chemical bonding for pair interactions with the COOP criterion [19] as implemented within the ASW method [START_REF] Methfessel | [END_REF]. They are based on non magnetic calculations in as far as the spin polarized calculation have the role of rigidly shifting the majority spins PDOS to lower energy and the minority spin PDOS to higher energy (Fig. 4b). The major part of the valence band (below E F ) is of bonding character with small antibonding COOP (negative magnitudes) from -4 eV up to E F arising from C1-C2 bonds due to their involving with U and Co bonding. The energy window is reduced with respect to the DOS panel in order to better show the pair interactions and because no significant bonding could be identified below -11 eV.

The major pair-bonding is in the energy range -8,-5 eV i.e. where C-p states are dominant. As expected the carbon pair (C1-C2) bonding prevails, followed by Co-C1 and U-C2 while the other interactions are of weaker intensity. U-Co bonding is also present especially in the energy range -4, -1 eV where Co-d states mix with uranium itinerant states below E F but it is weaker than the other intensities due to the large U-Co distance: ~3.6 Å. These bonding intensity follow from the crystal environments and the subsequent distances (Table 1). It is likely that the brittleness of the compound arises from such strong directional and site specific interactions.

Conclusions

In this work further evidence of the effects of C…C pairs in the building of the ternary carbide UCoC 2 was provided. The structural anisotropy arising from their alignment along the tetragonal c-axis leads to assign a high linear incompressibility as inferred from the energy-volume curves fitted with Birch EOS as well as from the set of elastic constants leading to assign a brittle like behavior. The chemical bonding is specific to carbon sites (C1 and C2) with cobalt and uranium respectively, and large C1-C2 and weak U-Co bonding. The nature of the C-C bond is found from the calculations of the frequencies to be closer to aliphatic C-C single bond. Further it is found that a dipole moment is developed due to the differentiated bonding of C1 and C2 respectively with Co and U. From energy differences the ground state is predicted ferromagnetic with major contribution to the magnetization arising from uranium spin polarization.

Future works will be devoted to experimental characterizations of the magnetic properties aided by full relativistic calculations. 1.

C 11 =

 11 C 22 = 358; C 12 = 241; C 13 = 149; C 33 = 530; C 44 = 206 and C 66 = 144. All C ij are positive and their combinations: C 11 > C 12 , C 11 C 33 > C 13 and (C 11 +C 12 )C 33 > 2C 13 2 , obey the rules pertaining to the mechanical stability of the chemical system. C 33 is much larger than C 11 = C 22 , in agreement with the high linear incompressibility along the c axis, i.e. in the direction of C…C pairs.The bulk (B V ) and shear (G V ) modules following Voigt are formulated as:B V = 1/9 {2(C 11 + C 12 ) + 4C 13 + C 33 } G V = 1/30 {12C 44 + 12C 66 + C 11 + C 12 + 2C 33 -4C 13 } The resulting numerical values are: B V = 258 GPa and Gv = 176 GPa. It is interesting to see that the bulk modulus obtained with the formula based on the elastic constants comes close to the one obtained by Birch EOS. This validates a posteriori the two approaches in analyzing the mechanical properties of the compound.

  results of computed charge changes (Q) are such that they lead to neutrality when the respective multiplicities are accounted for: U: +1.78; Co: +0.40; C1: -0.95 and C2: -1.22. These original results show different behaviors of carbon at the C1 and C2 sites with resulting charge difference of  = |0.27|. A C…C bond dipole moment (in Debye unit) is obtained:  = d  = 1.47  0.27 = 0.4 D. Despite its small magnitude compared with ionic compounds such as KBr (K + Br -with  = 10 D), this stresses the role of differentiated interactions and selective charge transfers of C1 and C2 with Co and U respectively through the chemical bond.

Fig. 1 .

 1 Fig. 1. (color) UCoC 2 : Sketch of the crystal structure including two formula units showing the C…C (C1-C2) pairs aligned along the tetragonal c-axis. Atom labels are as in Table1.

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .

 2345 Fig. 2. UCoC 2 : Energy volume curves for isotropic and anisotropic volume compressions and corresponding fit parameters (insert) from Birch 3 rd order EOS.

Table 1

 1 Experimental and calculated structural parameters of UCoC 2

	UCoC 2	Exp. [7]	Calc. PAW-GGA	GGA+U (U= 4eV)
	P4/nmm (Orig. 1)
	a (Å)	3.497	3.473	3.481
	c (Å)	7.313	7.275	7.372
	V (Å 3 )	89.43	87.75	89.34
	Co (2a) 0 0 0
	U,C1,C2 (2c) 0 ½ z
	z(U)	0.321	0.319	0.321
	z(C1)	0.8416	0.843	0.849
	z(C2)	0.6396	0.640	0.649
	Short distances (Å)
	U-C2	2.33	2.34	2.41
	Co-C1	2.09	2.08	2.07
	C1-C2	1.48	1.47	1.48
	E (eV) / 2 FU	-	-74.97	-67.47
	FU: formula unit.

Acknowledgements: Parts of the calculations were carried out on AVAKAS cluster of MCIA-Mesocenter of the University Bordeaux 1.