Tryfon Antonakakis

Fadi Baïda

Abderrahmane Belkhir

Kirill Cherednichenko

Shane Cooper

Richard Craster

Guillaume Demesy

John Desanto

Gérard Granet
email: gerard.granet@lasmea.univ-bpclermont.fr

Editor Evgeny Popov

Boris Gralak

Sébastien Guenneau

Daniel Maystre

André Nicolet

Brian Stout

Fréderic Zolla

Benjamin Vial

GratinGs: theory and numeric applications

Introduction

The C-method was born in the eighties in Clermont-Ferrand , France, from the need to solve rigorously diffraction problems at corrugated periodic surfaces in the resonance regime [START_REF] Chandezon | A new theoretical method for diffraction gratings and its numerical application[END_REF], [START_REF] Chandezon | Multicoated gratings: a differential formalism applicable in the entire optical region[END_REF], [START_REF] Li | Rigorous and efficient grating-analysis method made easy for optical engineers[END_REF]. The main difficulty of such problems is the matching of boundaries conditions. It is obvious that any method aimed at solving Maxwell's equation is all the more efficient since it is able to fit the geometry of the problem. For that purpose, Chandezon et al introduced the so called translation coordinate system deduced from the Cartesian coordinate system x, y, z) by the relations x = x 1 , y = x 2 , z = x 3 + a(x 1 ) where a(x 1 ) is a continuously differentiable function describing the surface profile. Hence since the boundary of the physical problem coincides with coordinate surfaces, writing boundary conditions is as simple as it is for classical problems in Cartesian, cylindrical, or spherical coordinates . This is the first ingredient of C-method. The second one is to write Maxwell's equation under the covariant form. This formulation comes from relativity where the use of curvilinear non orthogonal coordinate system is essential and natural. The main feature of this formalism is that Maxwell's equations remain invariant in any coordinate system, the geometry being shifted into the constitutive relations. Chandezon et al derived their 3D formulation from the general 4D relativistic Post's formalism [START_REF] Post | Formal structure of electromagnetics[END_REF] and evidently used tensorial calculus. Although it is with no doubt the most elegant and efficient way to deal with electromagnetic in general curvilinear coordinates it is also probably the reason why the theory appeared difficult to understand to many scientists. The third ingredient of C-method is that it is a modal method. This nice property is linked with the translation coordinate system in which a diffraction problem may be expressed as an eigenvalue eigenvector problem with periodic boundary conditions. The last feature of C-method is the numerical method of solution. The matrix operator is obtained by expanding field components into Floquet-Fourier harmonics and by projecting Maxwell's equations onto periodic exponential functions. The above four features may be resumed by saying that C-method is a curvilinear coordinate modal method by Fourier expansion [START_REF] Edee | Chandezon On the Rayleigh-Fourier method and the Chandezon method: Comparative study[END_REF]. Since the original papers, The C-method has gone through many stages of extension and improvement. The original theory was formulated for uncoated perfectly conducting gratings in classical mount. Various authors extended the method to conical diffraction mountings [START_REF] Plumey | Differential covariant formalism for solving Maxwell's equations in curvilinear coordinates:oblique scattering from lossy periodic surfaces[END_REF], [START_REF] Li | Multilayer-coated diffraction gratings: differential method of Chandezon et al. revisited[END_REF]. Granet et al [START_REF] Granet | Scattering by a periodically corrugated dielectric layer with non identical faces[END_REF], Li et al [START_REF] Li | some topics in extending the C method to multilayer-coated gratings of different profiles Pure[END_REF] and Preist et al [13] allowed the various profiles of a stack of gratings to be different from each other, although keeping the periodicity. Solving the vertical faces case in a simple manner, Plumey et al [START_REF] Plumey | Coordinate transformation method as applied to asymmetric gratings with vertical facets[END_REF] have shown that the method can be applied to overhanging gratings Preist et al obtained the same results by applying the usual coordinate transformation to oblique coordinates [START_REF] Preist | Optical response of blazed and overhanging gratings using Chandezon transformations[END_REF]. In the numerical context, Li [START_REF] Li | Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings[END_REF] and Cotter et al [START_REF] Cotter | scattering matrix approach to multilayer diffraction[END_REF] improved the numerical stability of the C-method by using the S-matrix propagation algorithm for multilayer gratings. It is seen that C-method has been applied to a large class of surface relief gratings and multilayer coated gratings. The key point of C-method is the joint use of curvilinear coordinates and covariant formulation of Maxwell's equations. All the new developments in the modelling of gratings like Adaptive Spatial Resolution [START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF], [START_REF] Granet | Reformulation of the coordinate transformation method through the concept of adaptive spatial resolution. Application to trapezoidal gratings[END_REF], [START_REF] Granet | Analysis of strip gratings using a parametric modal method by Fourier expansions[END_REF], and Matched coordinate [START_REF] Th | Matched coordinates and adaptive spatial resolution in the Fourier modal method[END_REF]derive from this fundamental observation.

C-Method

In Euclidean space with origin O and basis vector e x , e y , , e z , , let us consider an infinite cylindrical surface (Σ) whose elements are parallel to the y axis. This surface separates two linear homogeneous and isotropic media denoted (1) and [START_REF] Chandezon | Multicoated gratings: a differential formalism applicable in the entire optical region[END_REF]. In Cartesian coordinates such a surface can be described by equation z = a(x). Any electromagnetic field interacting with this particular geometry satisfies some boundary conditions. For instance, the tangential components of the electric field vector and the normal component of the displacement field vector are continuous at the surface. The point is that boundary conditions involve quantities that obviously depend on the position at which they are considered on the surface. We are thus led to look for a coordinate system which fits the problem and makes it more readily solvable than it is in a Cartesian framework. The so-called translation coordinate system (x 1 , x 2 , x 3 ) introduced by Chandezon and defined from the Cartesian coordinate system by the direct transformation (curvilinear coordinates to Cartesian coordinates) :

x = x 1 , y = x 2 , z = x 3 + a(x 1 ) (8.1)
or the inverse transformation (Cartesian coordinates to curvilinear coordinates):

x 1 = x, x 2 = y, x 3 = z -a(x) (8.2)
is one such system. It makes the surface (Σ) coincide with the coordinate surface x 3 = 0. A point M(x, y, z = a(x) at the surface (Σ) is now referenced by the triplet (x 1 , x 2 , 0). The coordinate surface x 3 = x 3 0 is obtained by translating each point at surface (Σ) with vector x 3 0 e z , hence the name given by Chandezon to this particular coordinate system: translation coordinate system. The change of coordinates may also be considered as a change of variable. This view . point allows a better understanding of the numerical behaviour of C-method and its connection with Rayleigh expansions. There is actually no difference in the way of deriving the elementary solutions of the scalar Helmholtz equation in Cartesian coordinates or in translation coordinate systems. Both are eigenvectors of an eigenvalue problem with pseudo-periodic boundary conditions. In both cases, the operator eigenvalue problem is transformed into a matrix eigenvalue problem thanks to the Galerkin method with pseudo periodic functions as expansion and test functions. Hence solving the scalar Helmholtz equation in any coordinate system is the very first step when implementing C-method. In the next paragraphs we shall focus on this issue before solving a grating problem.

Modal equation in the Cartesian coordinate system

Consider an homogeneous region with relative, possibly complex, permittivity, ε. In the harmonic regime with a time dependence of exp(-iωt), it is possible to construct general solutions to the field equations once we have general solutions to the scalar Helmholtz equation. So as a first task, we are going to investigate elementary solutions to the Helmholtz equation written in the translation coordinate system. Let us start from Cartesian coordinates in which 2D scalar Helmholtz equation is

∂ 2 x + ∂ 2 z + k 2 F = 0 (8.3)
where k = ω √ µ 0 ε is the wave-number. The coefficients of the Helmholtz equation are independent of z so we seek solutions of the form F (x, z) = exp(iγz)F(x) The Helmholtz equation becomes:

(∂ 2 x + k 2 )F(x) = γ 2 F(x) (8.4)
Function F(x) is thus an eigenmode of equation (8.4). The requirement that the eigenmodes satisfy the pseudo-periodicity condition F(x

+ d 1 ) = exp(iα O d 1 )F(x)
is automatically fulfilled by their expansion into Floquet-Fourier series:

F(x) = +∞ ∑ -∞ F m exp(iα m x) (8.5) α m = α 0 + mK 1 , K 1 = 2π d 1 , m ∈ N
and α 0 is some real parameter. By introducing (8.5) into (8.3) and by projecting onto pseudo-periodic functions exp(iα n x), one obtains the matrix equation:

γ 2 F = k 2 I -α F (8.6)
where F is a column vector whose elements are the F m and α is a diagonal matrix w hose elements are the α m and I is the identity matrix. The solution to the above matrix eigenvalue equation is of course trivial since the matrix is diagonal. Let us introduce subscript q to number the eigenvalues and the eigenfunctions. The eigenvalues γ q are deduced from their squared number:

γ 2 q = k 2 -α 2 q (8.7)
and the eigenvectors are determined by F mq = δ mq where δ mq is the Kronecker symbol. The square root of γ 2 q is defined as follows:

γ q =        γ 2 q if γ 2 q ∈ R + -γ 2 q if γ 2 q ∈ R - γ 2 q 1/2 with positive imaginary part if γ 2 q ∈ C (8.8)
Finally , F (x, z) can be represented by superposition of eigenmodes

F (x, z) = F + (x, z) + F -(x, z) (8.9) with F + (x, z) = q=+∞ ∑ q=-∞ A + q exp(iγ q z) m=+∞ ∑ m=-∞ δ mq exp(iα m x) (8.10) F -(x, z) = q=+∞ ∑ q=-∞ A - q exp(-iγ q z) m=+∞ ∑ m=-∞ δ mq exp(iα m x) (8.11)
There are two sets of modes, the number of which are equal: those propagating or decaying in the positive direction of z and those propagating or decaying in the opposite direction.

We denote these modes by superscript + and -respectively. The z dependence of an eigenmode is determined by function exp(iγ p z) By increasing z to z + ∆z , exp(iγ q z) is multiplied by exp(iγ q ∆z) = exp(iℜ(γ q )∆z) × exp(-ℑ(γ q )∆z) The real eigenvalues have ℑ(γ q ) = 0 and correspond therefore to forward modes if ℜ(γ q ) > 0 or backward modes if ℜ(γ q ) < 0. The complex eigenvalues modes have a non-zero imaginary part and possibly also a non-zero real part. The associated eigenmodes decay forward if ℑ(γ q ) > 0 or backward if ℑ(γ q ) > 0. These expansions are known as Rayleigh expansions; they are linear combination of eigenvectors that we call hereafter Rayleigh eigenvectors R q :

R q (x) = m=+∞ ∑ m=-∞ δ mq exp(iα m x) (8.12)
In Cartesian coordinates, the solutions to the Helmholtz equations may be regarded as the eigenvectors of a matrix equation. The eigenvalues γ m are determined by the periodic lateral boundary conditions of the problem and are obtained analytically since the matrix is diagonal. The translation coordinate system preserves the z translation symmetry and also periodic lateral boundary conditions. We may then expect a great formal similitude between solutions obtained in each coordinate system.

Modal equation in terms of the new variables

In this section, we derive the master equation of C-method by considering the change of coordinates as a change of variables. For the change of variables x 1 = x, x 2 = y, x 3 = za(x) the chain rule for derivatives has the form:

   ∂ x = ∂ 1 -ȧ∂ 3 ∂ y = ∂ 2 ∂ z = ∂ 3 (8.13)
Substituting the derivatives (8.13) into (8.3) gives:

(1 + ȧ ȧ)∂ 2 3 -ȧ∂ 1 ∂ 3 -∂ 1 ȧ∂ 3 + ∂ 2 1 + k 2 F (x 1 , x 3 ) = 0 (8.14)
the solution of which are the same as the solutions of (8.13) expressed in terms of the new variables.

F + a (x 1 , x 3 ) = F + (x = x 1 , x 3 = z -a(x)) = q=+∞ ∑ q=-∞ A + q exp(iγ q x 3 ) m=+∞ ∑ m=-∞
δ mq exp(iγ q a(x 1 )) exp(iα m x 1 ) (8.15)

F - a (x 1 , x 3 ) = F -(x = x 2 , x 3 = z -a(x)) = q=+∞ ∑ q=-∞ A - q exp(-iγ q x 3 ) m=+∞ ∑ m=-∞
δ mq exp(-iγ q a(x 1 )) exp(iα m x 1 ) (8.16)

The subscript a indicates the profile dependence of function F . We call function exp ±iγ q a(x 1 ) exp iα q x 1 the generalized Rayleigh eigenvector of order q. It is nothing more than plane wave exp(±iγ q z)exp(iα q x 1 ) expressed in terms of the new variables x 1 and x 3 and is closely linked with function a(x 1 ). Let us denote it R ± a,q : R ± a,q = exp ±iγ q a(x 1 ) exp(iα

q x 1 ) = m=+∞ ∑ m=-∞ R a mq exp(iα m x 1 ) (8.17)
It is assumed so far that a(x 1 ) is periodic with period d 1 hence:

exp(±iγ q a(x 1 )) = p+∞ ∑ p=-∞ L ± p exp i2π px 1 d 1 (8.18)
with:

L ± p = 1 d 1 d 1 0 exp ±iγ q a(x 1 ) exp -i2π px 1 d 1 dx 1 (8.19)
In physical space, the generalized Rayleigh eigenvectors result from the product of a periodic function with a pseudo-periodic one. Thus, in Fourier space, the spectrum of the qth generalized Rayleigh eigenvector is obtained by translating the spectrum of function exp(±iγ q z) with vector 2πq/d 1 , that is: R ± a mq = L ± a m-q (8.20)

Finally:

F + a (x 1 , x 3 ) = q=+∞ ∑ q=-∞ A + q exp(iγ q x 3 ) m=+∞ ∑ m=-∞ L + a m-q exp(iα m x 1 ) (8.21) F - a (x 1 , x 3 ) = q=+∞ ∑ q=-∞ A - q exp(-iγ q x 3 ) m=+∞ ∑ m=-∞ L - a m-q exp(iα m x 1 ) (8.22)
Functions (8.21) and (8.22) give the general solution to (8.14). Indeed each element of this solution is a generalized Rayleigh eigenvector associated to index q such that γ 2 q + α 2 q = k 2 and thus satisfies (8.14). The reason for that is obvious. It is obtain from (8.12) in which we have introduced the same change of variable as the one that has allowed us to get (8.14) from (8.3). From a practical view point, one can only manipulate finite size expansions and it does not make sense to speak of R ± a,q (x 1 ). That is why one may wonder if a generalized Rayleigh eigenvector is still a valid a solution of (8.14) when only a finite number of spatial Fourier harmonics is retained to represent it. Let us assume for a while the answer is yes and examine the involvements of such a claim. Introducing an integer M, hereafter denoted truncation number, and letting m run from -M to M the truncated generalized Rayleigh eigenvector writes:

R ±(M) a,q (x 1 ) = m=+M ∑ m=-M R ± a,mq exp(iα m x 1 ) (8.23)
Substituting ∂ 3 with iγ q , we have:

-(1 + ȧ ȧ)γ q Ṙ±(M) a -i( ȧ∂ 1 + ∂ 1 ȧ)γ q R ±,(M) a + (∂ 2 1 + k 2 )R ±,(M) a = 0 (8.24)
where Ṙ±(M) a denotes γ q R ±,(M) a

. Replacing ȧ by the coefficients of its Fourier series ȧp and denoting ȧ the toeplitz matrix whose elements ȧmp are the ȧm-p , it is easy to see that the matrix form of relation 8.24 is:

k 2 I -α 2 0 0 I R ± aq Ṙ± aq = γ q -ȧα -α ȧ I + ȧ ȧ I 0 R ± aq Ṙ± aq (8.25)
where R ± aq and Ṙ± aq are column vectors formed by the 2M + 1 Fourier coefficients of R ±(M) aq Ṙ±(M) aq respectively. (8.25) shows that γ q and R ± aq Ṙ± aq are an eigenvalue and an eigenvector of the generalized matrix eigenequation Aψ = ρBψ. Since R ± a,q is an exact eigenvector of (8.14) its truncated part can only approximate the solution of (8.14) and consequently, mathematically speaking γ q cannot be an eigenvalue of (8.25). It follows that our claim was false. Nevertheless, elementary pseudo-periodic solutions to (8.14) do exist and we will derive them in the next paragraph.

Fourier expansion of elementary waves in the translation coordinate system

In this paragraph, we derive the generalized eigenvalue eigenvector matrix equation starting from (8.14) the only assumption being the pseudo periodicity of the field and we discuss the obtained solutions. First, the propagation equation is rewritten as a pair of first-order equations:

k 2 + ∂ 2 1 0 0 1 F ∂ 3 F = -∂ 3 ȧ∂ 1 + ∂ 1 ȧ 1 + ȧ ȧ -1 0 F ∂ 3 F (8.26)
The coefficients of this equation do not depend on x 3 which allows to write the x 3 dependence as exp(iρx 3 ). The parameter ρ depends on the boundary conditions that F(x 1 , x 3 ) has to satisfy along x 1 direction. For gratings, periodic with period

d 1 along x 1 , F (x 1 + d 1 , x 3 ) = exp(iα 0 d 1 )F (x 1 , x 3
) where α 0 is some real parameter. ∂ 3 F verifies of course the same prop- erty. The above requirements on the solution are all fulfilled by expanding function F and ∂ 3 F under the form:

F (x 1 , x 3 ) = exp(iρx 3 )F a (x 1 ) = exp(iρx 3 ) m=M ∑ m=-M F a,m exp(iα m x 1 ) (8.27) ∂ 3 F (x 1 , x 3 ) = exp(iρx 3 ) Ḟa (x 1 ) = exp(iρx 3 ) m=M ∑ m=-M Ḟa,m exp(iα m x 1 ) (8.28)
Introducing the above expansions into (8.26) and projecting the latter onto exp i2πnx 1 d 1 basis, we get the sought algebraic matrix eigenvalue equation from which eigenvalues ρ q and eigenvectors F a,q are readily obtained thanks to standard computer libraries:

k 2 I -α 2 0 0 I F a,q Ḟa,q = ρ a,q -ȧα -α ȧ I + ȧ ȧ I 0 F a,q Ḟa,q (8.29)
As in the Cartesian coordinate system, it is observed numerically that there are two sets of modes, the number of which are equal: those propagating or decaying in the positive x 3 direction and those propagating or decaying in the opposite direction. Furthermore, it has been shown numerically and analytically [8], that, as the truncation number increases, the computed real eigenvalues converge to the real Rayleigh eigenvalues ±γ q .

lim M→∞ ±ρ M a,q = ±γ R q (8.30)

In the above relation, we have added an extra subscript M to indicate the truncation dependence. Indeed, the truncation order M has to be chosen large enough so that the computed real eigenvectors coincide with a great accuracy with their Rayleigh counterpart. In that case, provided that the eigenvalues are not degenerated, up to a multiplicative constant coefficient, the associated computed eigenvectors tend to the corresponding plane waves expressed in terms of the new variables

(x 1 , x 2 , x 3 ). lim M→∞ F ±(M) a,q = R ± a,q (8.31)
Thus in the translation coordinate system defined by x 1 = x, x 2 = y, x 3 = za(x) as in the Cartesian coordinate system Oxyz, linear combinations of elementary solutions to the Helmholtz equation allow us to express electromagnetic field while giving it a physical meaning in terms of forward and backward waves.We write numerically the solution to the Helmholtz equation as:

F + a (x 1 , x 3 ) = ∑ q∈U + A + q exp(iρ + a,q x 3 )R +(M) a,q (x 1 ) + ∑ q∈V + A + q exp(iρ + a,q (x 3 )F + a,q (x 1 ) (8.32) F - a (x 1 , x 3 ) = ∑ q∈U + A - q exp(iρ - a,q (x 3 )R -(M) a,q (x 1 ) + ∑ q∈V + A - q exp(iρ - a,q (x 3 )F - a,q (x 1 ) (8.33)
with:

F ± a,q (x 1 ) = m=+M ∑ m=-M F ± a,mq exp(iα m x 1 ) (8.34)
U ± , V ± denote the sets of indices for the propagating and decaying orders in the positive and negative direction respectively.

U + = q/ℜ(ρ a,q ) > 0 and ℑ(ρ a,q ) = 0 U -= q/ℜ(ρ a,q ) < 0 and ℑ(ρ a,q ) = 0 (8.35)

V + = q/ℑ(ρ a,q ) > 0 V -= q/ℑ(ρ a,q ) < 0 (8.36)

Application to a grating problem

Let's come back to the one-dimensional grating problem. Consider the electromagnetic problem in which two homogeneous non magnetic media are separated by a cylindrical periodic surface with period d 1 which is invariant along the y axis in the Cartesian coordinate system Oxyz. Such a surface, described by equation z = a(x) is illuminated from above by a unit amplitude linear polarized monochromatic plane wave with vacuum wavelength λ 0 , angular frequency ω and vacuum wave number k 0 = 2π/λ 0 . The wave vector is inclined at θ to the Oz axis. Medium (1) and medium (2) have relative permittivity ε 1 and ε 2 respectively. Time dependence is expressed by the factor exp(-iωt) Such a problem is reduced to the study of the two fundamental cases of polarisation and the unknown function F (x, z) is the y component of the electric or the magnetic field for TE and TM polarization respectively. We solved half the problem since we already determined the general solution to the scalar Helmholtz equation as a linear combination of elementary waves the coefficients of which remain to calculate. The situation is very common in electromagnetic theory: the fields on both side of the grating are expanded in terms of the modes in the respective regions with unknown coefficients. A method of solution known as mode-matching method was developed in the context of guided waves in the micro-wave range.

The grating may be considered as a generalized multi-port whose inputs are excited by waves that propagate or decay towards it giving rise to a response at the outputs that consists of the waves that propagate or decay away from it [START_REF]Numerical Techniques for microwave and millimeter-wave passive structures[END_REF], [START_REF] Redheffer | On the relation of transmission line theory to scattering and transfer[END_REF]. The mode coupling is caused by the modulation of the interface and by the different constitutive parameters in either side of it. The so-called scattering matrix S a defined as

A (1)+ A (2)-= S a A (1)- A (2)+ (8.37)
provides a linear relation between the output and input coefficients. In a grating problem, the vector formed by the amplitudes of the incoming waves has only one non null component: that corresponding to the incident wave which was assumed enforced to one. The subscript a indicates that the S matrix depends on the profile function a(x). We call S a matrix an interface scattering matrix. The S a matrix is derived from boundary conditions at the surface x 3 = x 3 0 . The change of variable makes it easy to write them. We have solved the scalar Helmholtz equation, the scalar field being a field component tangent to the surface;indeed F coincides with H y and E y in TM polarisation and TE polarisation respectively. For simplicity, let us consider TE polarisation where the non null components of the electromagnetic field are E y , H x , H z . Boundary conditions require matching the tangential components of the magnetic and electric field. We have already derived one of them, E y , we have to derive the tangential component H t of the electric field given by :

H t = H.t (8.38)
where t is the unit vector at point P which is tangential to the grating profile function. It is defined in terms of the e x and e z Cartesian unit vectors by:

t = 1 √ 1 + ȧ2 (e x + ȧe z ) (8.39)
The square root in the denominator represents a normalization factor that can be omitted since at a given point, it is identical on both sides of the boundary surface. let us introduce G such that:

G = iZH.t (8.40)
Where Z = µ 0 /ε is the wave impedance. From Maxwells equation we have iωµ 0 H x = -∂ z E y and iωµ 0 H z = ∂ x E y thus :

G (x, z) = - 1 k (∂ z F (x, z) -ȧ∂ x F (x, z)) (8.41)
substituting ∂ 3 for ∂ z and ∂ 1 -ȧ∂ 3 for ∂ x we get:

G (x 1 , x 3 ) = - 1 k ((1 + ȧ ȧ)∂ 3 -ȧ∂ 1 ) F (x 1 , x 3 ) (8.42)
Similarly to F , G depends on x 3 as exp(iρx 3 ) and we may write:

G (x 1 , x 3 ) = exp(iρx 3 )G(x 1 ) (8.43)
We are now familiar with the operational rules that allow to associate in Fourier space a matrix with an operator. We have:

1 + ȧ ȧ → I + ȧ ȧ, ȧ∂ 1 → i ȧα (8.44)
From which we deduce:

ikG ± a = (I + ȧ ȧ) F ± a ρ a -ȧαF ± a (8.45)
where ρ is a diagonal matrix whose elements are the eigenvalues ρ a,q . Writing the continuity of F (1) and F (2) and G (1) /Z (1) and G (2) /Z (2) at x 3 = x 3 0 is straightforward and leads to the following expression of the S a matrix:

S a = F (1)+ a -F (2)- a G (1)+ a -G (2)- a -1 -F (1)- a F (2)+ a -G (1)- a G (2)+ a (8.46)
The knowledge of S a matrix allows to calculate the constant coefficients of outgoing waves. Since the spectrum of the solutions of the transformed Helmholtz equation include the generalized Rayleigh eigenvectors associated to real Rayleigh eigenvalues the efficiencies may be calculated in the very same way as in the Cartesian coordinate system.

R q = |A (1)+ q | 2 γ (1) q γ (1) 0 T p = |A (2)- p | 2 γ (2) p γ (1) 0 (8.47) with: γ (1) q = k 2 0 ε 1 -k 0 √ ε 1 sin θ + q 2π d 1 2 γ (2) p = k 2 0 ε 2 -k 0 √ ε 1 sin θ + p 2π d 1 2 (8.48)
The values of integers p and q are such that γ

(1)

q and γ (2)
p are real.

Implementation of C-Method

The main interest to consider C-Method through a simple change of variable is to make us understand its numerical link with Rayleigh expansions and to calculate efficiencies as in the Cartesian Coordinate system. Within the framework of translation coordinate systems, starting from Maxwell's equations written under the covariant form, we have shown that all tangential components of the field at surface S could be generated from the longitudinal covariant components along the axis of invariance. Moreover,these components are solutions of the scalar Helmholtz equation. Therefore, one clearly understands that finding the elementary solutions of the scalar Helmholtz equation is the kernel of C-method. To summarize we may enunciate the different steps for solving a grating problem with C-method:

• Define a translation coordinate system.

• Find the elementary waves of the Helmholtz equation. For that purpose use the Galerkin method with exp(iα n x 1 ) as expansion and test functions. Substitute the generalized Rayleigh eigenvectors for the computed eigenvectors associated to real eigenvalues. Sort the elementary waves into forward and backward waves.

• Write boundary conditions at surface Σ and calculate efficiencies as in the Cartesian coordinate system.

Various formulations of C-method

So far, the Helmholtz equation in the translation coordinate system was derived by using the chain rule for derivatives in the Helmholtz equation written in the translation coordinate system. In this section,we start from the covariant Maxwell's equations and we show that they lead to several operators one of them being the propagation equation. In a homogeneous isotropic medium with permittivity ε and permeability µ, with a time dependence exp(-iωt),the symmetrized Maxwell equations write:

ξ αβ γ ∂ β F γ = k √ gg αβ G α ξ αβ γ ∂ β G γ = k √ gg αβ F α (8.49)
where k = ω √ µε, the F γ and the G β are the complex amplitudes of the electric field and of a renormalized magnetic field respectively. We restrict our analysis to 1D problems in which both the geometry and the solution are independent of y. Practically this means that ∂ 2 is null as well as g 12 , g 21 , g 32 and g 23 . It follows that (8.49) decouple into two fully identical systems where the non null components are F 2 , G 1 , G 3 , and G 2 , F 1 , F 3 respectively.The first set of three components corresponds to T E polarisation, the second one to T M polarisation. Both polarisations obey the same first order differential equations system written hereafter for T E polarisation:

-∂ 3 F 2 = k √ gg 11 G 1 + √ gg 13 G 3 (8.50a) ∂ 1 F 2 = k √ gg 31 G 1 + √ gg 33 G 3 (8.50b) ∂ 3 G 1 -∂ 1 G 3 = k √ gg 22 F 2 (8.50c)
For T M polarisation, it is enough to permute F and G . Among the three components of each system, two play a particular role. Let us assume that x 3 = x 3 0 separates two isotropic homogeneous media. Then, in T E polarisation F 2 and G 1 / µ ε have to be continuous at surface

x 3 = x 3 0 .
The same conclusions holds for G 2 / µ ε and F 1 for T M polarisation. So, we have to solve (8.50) for the components labelled by two and by one. C-method is a Fourier based method which means that constitutive relations have to be written in Fourier space. In other words a matrix is to be associated to each element √ gg αβ of the constitutive tensors. The way for doing so should follow the so-called "Fourier factorization" rules derived by Li [START_REF] Li | Use of Fourier series in the analysis of discontinuous periodic structures[END_REF], [START_REF] Li | Oblique-coordinate-system-based Chandezon method for modeling one-dimensinally periodic, multilayer, inhomogeneous, anisotropic gratings[END_REF]. let us denote by ( √ gg αβ ) the matrix associated to coefficient √ gg αβ . According to Li's rules, the ( √ gg αβ ) write:

( √ gg 11 ) = 1 √ gg 11 -1 ( √ gg 13 ) = 1 √ gg 11 -1 g 13 g 11 ( √ gg 31 ) = g 31 g 11 1 √ gg 11 -1 ( √ gg 33 ) = 1 √ gg 11 + g 31 g 11 1 √ gg 11 -1 g 13 g 11 ( √ gg 22 ) = √ gg 22 (8.51)
The notation [ f ] designates the toeplitz matrix whose elements f mp are the f m-p elements of the Fourier series of function f (x 1 ). For the translation coordinate (x 1 , x 2 , x 2 ) such that x = x 1 , y = x 2 , z = x 3 + a(x 1 ), we have:

√ gg 11 → I √ gg 13 → -ȧ √ gg 31 → -ȧ √ gg 33 → [I + ȧ ȧ] √ gg 22 → I (8.52)
In Fourier space, the derivative operator ∂ 1 is associated to the diagonal matrix iα the elements of which are the iα m such that:

α m = α 0 + m 2π d 1 (8.53) Setting F 2 (x 1 , x 3 ) = m=+M ∑ m=-M F 2m (x 3 ) exp(iα m x 1 ) (8.54a) G 1 (x 1 , x 3 ) = m=+M ∑ m=-M G 1m (x 3 ) exp(iα m x 1 ) (8.54b) G 3 (x 1 , x 3 ) = m=+M ∑ m=-M G 3m (x 3 ) exp(iα m x 1 ) (8.54c) (8.54d)
We are now able to write (8.50) in Fourier space:

-∂ 3 F 2 = k √ gg 11 G 1 + √ gg 13 G 3 (8.55a) iαF 2 = k √ gg 31 G 1 + √ gg 33 G 3 (8.55b) ∂ 3 G 1 -iαG 3 = k √ gg 22 F 2 (8.55c)
where F 2 , G 1 G 3 are column vectors of size 2M + 1 whose components are the F 2m (x 3 ), G 1m (x 3 ), G 3m (x 3 ) respectively:

F 2 (x 3 ) = [F 2,-M (x 3 ), F 2,-M+1 (x 3 ), • • • , F 2,0 (x 3 ), • • • , F 2,M-1 (x 3 ), F 2,M (x 3 )] T (8.56a) G 1 (x 3 ) = [G 1,-M (x 3 ), G 1,-M+1 (x 3 ), • • • , G 1,0 (x 3 ), • • • , G 1,M-1 (x 3 ), G 1,M (x 3 )] T (8.56b) G 3 (x 3 ) = [G 3,-M (x 3 ), G 3,-M+1 (x 3 ), • • • G 3,0 (x 3 ), • • • , G 3,M-1 (x 3 ), G 3,M (x 3 )] T (8.56c)
where the exponent T is for the transposition.

Propagation equation in curvilinear coordinates

From Eqs (8.50a) and (8.50b), G 1 and G 3 may be expressed in terms of

F 2 kG 1 = - √ gg 33 ∂ 3 F 2 + √ gg 13 iαF 2 (8.57) kG 3 (x 3 ) = √ gg 31 ∂ 3 F 2 (x 3 ) + √ gg 11 iαF 2 (x 3 ) (8.58)
(8.50c), in which we substitute G 1 and G 3 with expressions (8.57) and (8.58), gives the propagation equation:

-α √ gg 11 α + ∂ 3 √ gg 33 ∂ 3 + iα √ gg 13 ∂ 3 + ∂ 3 √ gg 31 iα + k 2 √ gg 22 F 2 (
x 3 ) = 0 (8.59) which is rewritten as a pair of first-order differential equation as:

-i∂ 3 A F 2 (x 3 ) -i∂ 3 F 2 (x 3 ) = B F 2 (x 3 ) -i∂ 3 F 2 (x 3 ) (8.

60)

with:

A = α √ gg 13 + √ gg 13 α √ gg 33 I 0 (8.61) B = -α √ gg 11 α + k 2 √ gg 22 0 0 I (8.62)
Since the coefficients of matrices A and B are independent of variable x 3 , we may seek vectors F 2 (x 3 ), G 1 (x 3 ),G 3 (x 3 ) under the form:

F 2 (x 3 ) = F 2 exp(iρx 3 ) (8.63a) -i∂ 3 F 2 (x 3 ) = Ḟ2 exp(iρx 3 ) (8.63b) G 1 (x 3 ) = G 1 exp(iρx 3 ) (8.63c) G 3 (x 3 ) = G 3 exp(iρx 3 ) (8.63d)
This last step transforms (8.60) into a generalized eigenvalue eigenvector matrix equation:

Aρ F 2 Ḟ2 = B F 2 Ḟ2 (8.64)
It is then easy to check that (8.64) is the same as (8.29). After F 2 is determined, it remains to deduce G 1 from (8.57).

"Classical" C-method operator

We call "classical" operator the operator derived by Chandezon in his early work. From (8.50b) and taking into account (8.63) we find an expression for G 3 as follows:

G 3 = 1 k √ gg 33 -1 iαF 2 - √ gg 33 -1 √ gg 31 G 1 (8.65)
Substituting G 3 in (8.50a) and (8.50c) with the above expression yields:

   - √ gg 13 √ gg 33 -1 α ik √ gg 11 - √ gg 13 √ gg 33 -1 √ gg 31 -ik √ gg 22 - 1 k 2 α √ gg 33 -1 α -α √ gg 33 -1 √ gg 31    F 2 G 1 = ρ F 2 G 1 (8.66)

Multilayer grating

The extension of C-method to multilayer gratings is straightforward provided the interfaces which separate the layers share the same periodicity. It is just a generalization of the theory of planar stratified media. As a canonical case, let us consider a layer made of isotropic homogeneous media limited on the top by surface z = a j (x 1 ) and on the bottom by surface z = a j+1 (x 1 ) = a j (x 1 )t j (x 1 ). When t j (x 1 ) is constant the two surfaces are parallel to each other. In homogeneous media the field is a superposition of forward and backward waves. The only places where coupling occurs are the interfaces. Thus, we have to describe two different phenomena: on the one hand scattering at the interfaces and on the other hand propagation or attenuation in the layer. To summarize we assimilate an interface to a 4N-port local network (N = 4M + 1, M being the truncation number) and a layer to a multi-channel pipe connecting the 2N-ports of its input network and output network [START_REF] Redheffer | On the relation of transmission line theory to scattering and transfer[END_REF]. We have already defined interface scattering matrices which are local matrices in the sense they depend on the profile. In other words,in the context of C-method they depend on the coordinate system. Thus, for a layer bounded by two non parallel surfaces, we have to solve two eigenvalue problems for each surface which allows to calculate interface matrices S a j and S a j+1 . It remains to define and to calculate layer scattering matrices. Although two cases have to be considered according to whether the layer separates two identical surfaces or not, the line of reasoning is the same. As already mentioned, we have two coordinate systems such that z = x 3 j + a j (x 1 ) and z = x 3 j+1 + a j+1 (x 1 ). They are linked by the following relation:

x 3 j = x 3 j+1 + a j+1 (x 1
)a j (x 1 ) = x 3 j+1t j (x 1 ) (8.67) In medium j, located in between surfaces z = a j (x 1 ) and z = a j+1 (x 1 ), we may express the linear combination of forward and backard waves with coordinate x 3 j = 0 as local origin (that is z = a j (x 1 ) and write:

F ( j) a j (x 3 j , x 1 ) = ∑ q A ( j)+ j,q exp(iρ ( j)+ a j ,q x 3 j )F ( j)+
a j ,q (x 1 ) + ∑ q A ( j)j,q exp(iρ ( j)a j ,q x 3 j ))F ( j)a j ,q (x 1 ) (8.68)

In the same medium j we may also choose x 3 j+1 = 0 as local origin (that is z = a j+1 (x 1 )) which gives:

F ( j) a j+1 (x 3 j+1 , x 1 )) = ∑ q A ( j)+ j+1,q exp(iρ ( j)+ a j+1 ,q x 3 j+1 )F ( j)+ a j+1 ,q (x 1 )+ ∑ q A ( j)- j+1,q exp(iρ ( j)- a j+1 ,q (x 3 j+1 )F ( j)- a j+1 ,q (x 1 ) (8.
69) The layer is considered as a 4Nports which connects input waves F ( j)a j and F ( j)+ a j+1 to output waves F ( j)+ a j and F ( j)a j+1 , hence the definition of the layer S matrix:

  A ( j)+ j,q A ( j)- j+1,q   = S j, j+1   A ( j)- j,q A ( j)+ j+1,q   (8.70)
At the input of the layer, that is at x 3 j = 0, the outgoing waves correspond to the incoming wave of the output plane:

F ( j)+ a j (x 3 j = 0) = F ( j)+ a j+1 (x 3 j+1 = t j (x 1 )) (8.71)
Similarly, at the output of the layer, that is x 3 j+1 = 0, the outgoing waves correspond to the incoming waves of the input plane:

F ( j)- a j+1 (x 3 j+1 = 0) = F ( j)- a j (x 3 j = -t j (x 1 ) (8.72)
At this stage, we infer that layer S matrix looks like:

S j, j+1 = 0 P ( j)+ P ( j)- 0 (8.73)

The sought sub-matrices P ( j)+ , P ( j)-depend on whether the layer faces are parallel or not.

Layer with non parallel faces

Consider equation (8.75) and write it in terms of the eigenvectors of both coordinate systems:

∑ m ∑ q A ( j)+ j,q F ( j)+ a j ,mq exp(iα m x 1 ) = ∑ m ∑ q A ( j)+
j+1,q exp(iρ + a j +1,q t j (x 1 ))F

( j)+ a j+1,q exp(iα m x 1 ) (8.74)

The left hand side purely consists of a linear combination of eigenvectors expanded onto the exp(iα m x 1 ) basis whereas the right hand side consists of a linear combination of eigenvectors each of which being multiplied by a periodic functions of the x 1 variable. In order to get a matrix relation between the A ( j)+ j,q and the A ( j)+ j+1,q , we project (8.74) onto exp(iα m x 1 ). We get:

∑ m ∑ q A ( j)+ j,q F ( j)+ a j ,mq = ∑ m ∑ q A ( j)+ j+1,q F( j)+ a j+1,q exp(iα m x 1 ) (8.75)
with:

F( j)+ a j+1 m q = 1 d 1 d 1 0 ∑ l F ( j)+
a j+1 ,l q exp(iα l x 1 ) exp(iρ ( j)+ a j +1,q t j (x 1 )) exp(-iα m x 1 )dx 1 (8.76)

Then, the P ( j)+ matrix is readily obtained as

P ( j)+ = F ( j)+ a j -1 F ( j)+ a j+1 (8.77)
where

F ( j)+ a j
(respectively F ( j)+ a j+1 ) is the matrix formed by juxtaposition of vectors F ( j)+ a j ,q (respectively F( j)+ a j+1 q . Similarly we have:

∑ m ∑ q A ( j)- j+1,q F ( j)- a j+1 ,mq exp(iα m x 1 ) = ∑ m ∑ q A ( j)- j,q exp(iρ ( j)- a j ,q t j (x 1 ))F ( j)- a j,q exp(iα m x 1 ) (8.78) and ∑ m ∑ q A -( j) j+1,q F ( j)- a j+1 ,mq = ∑ m ∑ q A ( j)- j,q F( j)- a j,q (8.79)
with:

F( j)- a j m q = 1 d 1 d 1 0 ∑ l F ( j)-
a j ,l q exp(iα l x 1 ) exp(-iρ ( j)a j +1,q t j (x 1 )) exp(-iα m x 1 )dx 1 (8.80) from which we derive:

P ( j)-= F ( j)- a j+1
-1 F ( j)a j (8.81)

Layer with parallel faces

In that case, the two coordinate systems are identical and t j (x 1 ) is a constant. Equations (8.72) and (8.75) reduce to:

∑ m ∑ q A ( j)+ j,q F ( j)+ a j ,mq exp(iα m x 1 ) = ∑ m ∑ q A +( j)
j+1,q exp(iρ

( j)+ a j ,q t j )F ( j)+ a j ,mq exp(iα m x 1 ) (8.82) ∑ m ∑ q A ( j)- j+1,q F ( j)+ a j ,q = ∑ m ∑ q A
( j)j,q exp(-iρ ( j)a j ,q t j )F ( j)a j ,mq exp(iα m x 1 ) (8.83) from which we easily deduce:

A ( j)+ j,q = A ( j)+
j+1,q exp(iρ ( j)+ a j ,q t j ) or P ( j)+ = diag exp(iρ

( j)+ a j ,q t j ) (8.84) A ( j)- j+1,q = A ( j)- j,q exp(-iρ ( j)- a j ,q t j ) or P ( j)-= diag exp(-iρ ( j)- a j ,q t j ) (8.85)
It should be noted that when ρ ( j)+ a j ,q (respectively ρ ( j)a j ,q ) is complex valued, its imaginary part is negative (respectively positive) . Since t j is positive, exponential functions exp(±iρ ( j)± a j ,q t j ) associated to complex eigenvalues always decay when the layer thickness increases.

Combination of S matrices

The final step for analysing reflection and transmission by a layer is to combine the two interfaces S matrix and the layer S matrix. The tool for doing this is the Redheffer star product which gives the composition rules of two cascaded S matrices [26]. Consider two S matrices and partition them into four blocks: where U 1 = I -S 22 a j P ( j)+ S 11 a j+1 P ( j)--1 U 2 = I -S 11 a j+1 P ( j)-S 22 a j P ( j)+ -1

S 1 = S 11 1 S 12

Extensions of C Method

The key idea of C-method as applied to diffraction by surface-relief gratings is to map the surface of the grating to a plane. Until now, we have only described profiles under the form z = a(x). However, in Cartesian coordinates (x, y, z), a cylindrical surface whose generating line is parallel to the Oy axis may be described by the parametric equations:

x = f (x 1 ) z = g(x 1 ) (8.97)
where f and g are two continuous functions. Now consider the following relations:

x = f (x 1 ) + c 1 x 3 y = x 2 z = g(x 1 ) + x 3 (8.98)
where c 1 is a real constant. They define an additive change of coordinates whose metric tensor is given by:

g i j =   ∂ 2 1 f + ∂ 2 1 g 0 c 1 ∂ 1 f + ∂ 1 g 0 1 0 c 1 ∂ 1 f + ∂ 1 g 0 1 + ∂ 2 1 g   (8.99)
Actually, the above matrix corresponds to a change of coordinates provided the Jacobian determinant J of the transformation does not go to zero.

J = ∂ 1 x ∂ 3 x ∂ 1 z ∂ 3 z = ∂ 1 f c 1 ∂ 1 g 1 = ∂ 1 f -c 1 ∂ 1 g (8.100)
More over, the metric tensor is independent of coordinate x 3 which means there exists a translation symmetry along x 3 axis. Hence Equations(8.98) define in a general way translation coordinate systems which allow to solve new classes of problems.

Oblique transformations

In Cartesian coordinates, usual coordinates lines of a plane are two straight lines orthogonal to each other. One can also imagine having straight lines which make an angle different from π/2. Consider the straight line ∆ given by z = tan(φ )x and let us call φ the obliquity angle. The following sets of relations define a coordinate system (x 1 , x 3 ) in which lines parallel to ∆ are coordinate lines x 1 = constant and lines x 3 = constant remain parallel to Ox

x = x 1 + 1 tan φ x 3 z = x 3 (8.101)
Such oblique transformation allow to model an extended class of surface shapes which would otherwise be numerically inefficient (very blazed gratings) or even impossible like overhanging gratings. As an illustrative example, consider in the coordinate system (x 1 , x 3 ) the symmetric triangular function. Using an oblique transformation one gets:

t(x 1 ) =      2x 1 0 < x 1 < .5 2(1 -x 1 ) .5 < x 1 < 1 0 elsewhere (8.102)
x = x 1 + 1 tan φ (x 3 + t(x 1 )) z = x 3 + t(x 1 ) (8

Stretched coordinates

The essence of C-method is to choose a coordinate system that facilitate the solution of a given problem. Oblique transformations are a typical example of the usefulness of this technique. Indeed they provide an easy and elegant way to handle gratings with one vertical facet and also overhanging gratings. Similarly, we have believed for a long time that sharp edges were an intrinsic limitation of the C method. Actually, it turns out that transformations which stretch coordinates around the edges overcome the problem. With C-Method, the solution of Maxwell's equations is reduced to the solution of an algebraic eigenvalue problem in discrete Fourier space. The derivation of the matrix operator involves two steps: (1) The electromagnetic field is expanded into Floquet-Fourier series, and (2) the derivative of the grating profile function is expanded into Fourier series. When the latter function is discontinuous, the Fourier method is known to converge slowly.This weakness remains even when the correct Fourier factorization of products of discontinuous periodic functions, as given by Li, is applied. The reason for slow convergence is that the spatial resolution of the Fourier expansion remains uniform within a grating period whatever the grating profile function may be. On the contrary, stretched coordinates allow a mapping of space that increases spatial resolution around the discontinuities of the derivative of the profile function. For this reason the technique is known as adaptive spatial resolution.

Parametric C-method

Whether for mandatory reasons as is the case for overhanging gratings or simply to improve convergence speed, the most general representation of a one dimensional profile happens to be a parametric one. Adding an additional degree of freedom with an obliquity angle, a class of translation coordinate systems has the form given by (8.98). Due to the translational symmetry along vector e 3 = c 1 e x + e z , a numerical solution in terms of eigenvectors and eigenvalues is possible. Equations (8.98) describe a coordinate system where coordinates lines x 3 = constant coincide with functions which are periodic with period d 1 along direction Ox. Compared to the non-oblique coordinate system, the period d 1 and the direction of periodicity remain unaffected by the introduction of parameter c 1 . Thus, assuming an incident plane wave vector k such that k.e x = α 0 the x 1 dependence is of the form exp(iα m x 1 ) with α m = α 0 + m 2π d 1 . So we have all the ingredients to determine the matrix from which eigenvectors and eigenvalues will be sought. In Fourier space, the matrices associated to the elements of the metric tensor are:

√ gg 11 = (c 2 1 + 1) ḟ -c 1 ġ -1 √ gg 13 = ḟ -c 1 ġ -1 c 1 ḟ + ġ √ gg 31 = c 1 ḟ + ġ ḟ -c 1 ġ -1 √ gg 33 = ḟ -c 1 ġ + (c 2 1 + 1) c 1 ḟ + ġ ḟ -c 1 ġ -1 c 1 ḟ + ġ √ gg 22 = ḟ -c 1 ġ (8.104)
where ḟ and ġ designates the toeplitz matrices formed by the elements of the Fourier series of ∂ 1 f and ∂ 1 g respectively.

Plane waves and parametric C-method

More over since the physics remains the same compared to non-oblique translation coordinate systems, eigenvectors separate into forward and backward waves as was already the case:

F (x 1 , x 3 ) = ∑ q A ± q F ± q (x 1 )exp(iρ ± q x 3 ) (8.105)
As in the classical translation coordinate system, we substitute the computed propagative forward and backward eigenvectors with the corresponding transformed plane waves . Consider plane waves exp(iα n x) exp(±iγ n z) such that ±γ n ∈ R. Taking into account (8.98),their expression in oblique coordinates is:

exp i(α n c 1 ± γ n )x 3 exp i(α n f (x 1 ) ± γ n g(x 1 ) (8.106)
hence the following correspondences between propagative waves in Cartesian coordinates and their computed counterparts in oblique coordinates:

ρ ±,(M) (p),n ←→ (±γ n + α n c 1 ) ; F ± (p),n (x 1 ) ←→ exp iα n f (x 1 )±iγ n (g(x 1 ) (8.107)
We have added an extra subscript (p) and a superscript (M) to indicate that we only care about the above correspondence for propagative waves and that ρ depends on the truncation number.

Illustrative example

Consider a right angled triangular profile whose base is aligned on Ox. Other parameters are period d 1 and height h. It is illuminated by a plane wave inclined at θ to the Oz axis. In the context of C-method we ask ourselves which coordinate system choosing for modelling diffraction by such a grating. Here the main difficulty comes from the vertical facet located at x = d 1 . The operator associated with C-method involves the derivative of the profile function.

With a description of the profile by a function of the kind z = a(x), the derivative is constant and everything happens as if the vertical did not exist. Should the vertical be replaced by a very sloping facet, then a highly located and large discontinuity in the derivative would appear. None of the situation is satisfactory. An easy way to overcome the problem consists in introducing an oblique coordinate system in which the vertical is transformed into a straight line with a "reasonable" slope. Actually, doing so amounts to parametrizing the profile in the Cartesian coordinate system.

Obliquity angle and parametrization of the profile

Since one of the facets of the grating is vertical, an inclined coordinate system is needed. On the one hand, the parameter c 1 is linked to the obliquity angle φ by c 1 = 1/ tan φ and on the other hand, according to (8.100) it should satisfy the constraint 1c 1 ∂ x a > 0 . Hence, in principle φ may be any angle such that tan φ < h/d Let t 1 be tan(φ ). On the first facet we have:

x = x 1 + 1 t 1 y, y = h d x (8.108)
and on the second one

x = d, d = x 1 + 1 t 1 y (8.109)
Thus the parametrization of the profile is:

f (x 1 ) = d d -t 1 h x 1 g(x 1 ) = h d -t 1 h x 1 if x 1 x 1 0 f (x 1 ) = d g(x 1 ) = 1 t 1 (x 1 -d) if x 1 > x 1 0 (8.110)
with :

x 1 0 = d 1 - 1 t 1 h d (8.111)
Now that we have parametrized the profile, it remains to define a translation direction. The direction which served to parametrize the profile is a natural choice although not mandatory. Once more, the only constraint is that the tangent of the chosen obliquity angle is smaller than h/d.

Stretched coordinates and parametrization of the profile

At x 1 = 0 and x 1 = x 1 0 , the parametric functions f (x 1 ) and g(x 1 ) have jumps which can be reduced if one introduces an additional change of coordinates aimed at increasing spatial resolution around these points. Let x 1 be a function of a new variable u: x 1 = s(u). The chain rule for derivative gives :

∂ u x = ∂ 1 f (x 1 (s(u))∂ u s, ∂ u y = ∂ 1 g(x 1 (s(u))∂ u s (8.112)
Compared to the initial parametrization, spatial resolution is modulated by the multiplicative factor ∂ u s. The smaller the latter, the higher the spatial resolution. A possible stretching function is as follows :

s(u) =          u - ηx 0 2π sin 2πu x 0 if 0 u < x 1 0 (u -x 0 ) - η(d 1 -x 1 0 ) 2π sin 2π(u -x 0 ) d 1 -x 1 0 if x 1 0 u < d 1 (8.113)
The parameter η between zero and one controls the density of coordinate lines around the transition points. It allows to stretch space thinner where discontinuities of coefficients in Maxwell's equations occur. The larger η, the smaller ∂ u s and thus the higher the spatial resolution. In principle the parameter η does not have to reach one because, in that case, the Jacobian would be zero. Figure (8.6) shows four possible parametrization of the considered right angle triangle. Case (a) corresponds to the usual representation z = a(x). Finally, figure (8.7) shows the speed of convergence of the specular reflected order for a perfectly conducting right angle triangular profile for two different parametrizations. It has to be emphasized that modelling this kind of profile is out of reach for the "classical" C-method since it has a vertical facet. In Cartesian coordinates we deal with three mutually perpendicular families of planes: x=constant, y=constant, z=constant. Imagine that we superimpose on this system three other families of surfaces. We may reference any variable point M by the intersection of three planes in Cartesian coordinates , ie by the triplet (x, y, z) or as the intersection of the three surfaces that form our new, curvilinear coordinates. Describing the curvilinear coordinates surfaces by x 1 = constant,x 2 = constant,x 3 = constant we may identify our point by the triplet x, y, z as well as by x 1 , x 2 , x 3 . This means that in principle we may define a curvilinear coordinate system from the Cartesian system (x, y, z) by:

x = x 1 ′ = x 1 ′ (x 1 , x 2 , x 3 ), y = x 2 ′ = x 2 ′ (x 1 , x 2 , x 3 ), z = x3 = x 3 ′ (x 1 , x 2 , x 3 ) (8.114)
or by the inverse relations

x 1 = x 1 (x 1 ′ , x 2 ′ , x 3 ′ ), x 2 = x 2 (x 1 ′ , x 2 ′ , x 3 ′ ), x 3 = x 3 (x 1 ′ , x 2 ′ , x 3 ′ ) (8.115) x 1 ′ , x 2 ′ , x 3 ′ respectively x 1 , x 2 ,
x 3 are regarded as independent and continuously differentiable functions of x 1 , x 2 and x 3 respectively x 1 ′ ,x 2 ′ ,x 3 ′ . let M denote a variable point referenced by the rectangular coordinates(x, y, z).At M the so-called natural referential(M, e 1 , e 2 , e 3 )is defined by the the following basis vectors:

e α = β ′ =3 ∑ β ′ =1 ∂ x β ′ ∂ x α e β ′ (8.116)
with e 1 ′ = e x , e 2 ′ = e y , e 3 ′ = e z , e x , e y and e z being the unit vectors of an orthogonal Cartesian referential. In a similar way we may write

e α ′ = β =3 ∑ β =1 ∂ x β ∂ x α ′ e β (8.117) Moreover introducing Λ β ′ α = ∂ x β ′ ∂ x α
and Einsteins' summation convention Eq(8.116) and Eq(8.117) write: is the Jacobian matrix J of the change of coordinates. Since functions x 1 ′ , x 2 ′ ,x 3 ′ are independent J is inverible and its inverse is formed by the coefficients

e α = Λ β ′ α e β ′ e α = Λ β α e β ′ ( 8
Λ β α ′ J =   Λ 1 ′ 1 Λ 2 ′ 1 Λ 3 ′ 1 Λ 1 ′ 2 Λ 2 ′ 2 Λ 3 ′ 2 Λ 1 ′ 3 Λ 2 ′ 3 Λ 3 ′ 3   J -1 =   Λ 1 1 ′ Λ 2 1 ′ Λ 3 1 ′ Λ 1 2 ′ Λ 2 2 ′ Λ 3 2 ′ Λ 1 3 ′ Λ 2 3 ′ Λ 3 3 ′   (8.119)
One can also define basis vectors e α that are normal to coordinate surfaces x α = constant by

e α = ∂ x α ∂ x α ′ e α ′ , with e α ′ = e α ′ (8.120)
The vectors e α and e β form a set of reciprocal basis with e α .e β = δ β α , where δ β α is the Kronecker delta. The representation of any vector A in one of these bases is: The A α and the A α are the contravariant components and the covariant components of vector A respectively. The nullity of a component of vector A may be geometrically interpreted as follows:

A = A α e α = A α e α (8.121)
A α = 0: A is orthogonal to the tangent at point M to the coordinate line x α A α = 0: A belongs to the tangential plane at point M to coordinate surface x α In normalized orthogonal Cartesian coordinates differentiate contravariant and covariant components of a vector is generally not necessary. The Jacobian matrix allows to express the Cartesian components A α ′ = A α ′ of vector A in terms of its local contravariant components A α = A.e α or covariants components A α = A.e α :

A α ′ = A α ′ = Λ α ′ α A α or A α = Λ α α ′ A α ′ (8.122) A α ′ = Λ α α ′ A α or A α = Λ α ′ α A α ′ (8.123)
The quantities g αβ = e α .e β = Λ α ′ α Λ β ′ β g α ′ β ′ (8.124) define the metric of the coordinate system. In matrix form we have: We have seen that the natural referentiel gives the tools to easily maipulate tangential and normal components of a vector field.Therefore, writting boundary conditions at a surface should be straightforward. We need now to express Maxwell's equation in the new coordinate system.

g αβ = J t J (8.
For that purpose, we may follow a tensorial approach or stay at an elementary level and make a simple change of coordinates and components in the usual Maxwell's equations. We present briefly both points of view. A time dependence of the form exp(-iωt) is assumed.

Vectorial approch

Let us start from one of the Maxwell's curl equation written in the Cartesian coordinate system and in an homogeneous medium with permittivity ε and permeability µ ξ α ′ β ′ γ ′ ∂ β ′ H γ ′ = -iωεE α ′ (8.128)

where ξ α ′ β ′ γ ′ stands for the Levi-Civita indicator :

ξ α ′ β ′ γ ′ =    1 for α ′ β ′ γ ′ =
H γ ′ = Λ β γ ′ H β ξ α ′ β ′ γ ′ Λ α β ′ ∂ α Λ β γ ′ H β = -iωεE α ′ (8.130)
The left hand side of the above equation is equal to:

ξ α ′ β ′ γ ′ Λ α β ′ Λ β γ ′ ∂ α H β + ξ α ′ β ′ γ ′ Λ α β ′ ∂ α Λ β γ ′ H β (8.131)
on the one hand we have

Λ α β ′ ∂ α Λ β γ ′ = ∂ β ′ Λ β γ ′ (8.132)
and the other hand this term is symmetrical with respect to β ′ and γ ′ . Thus by applying the operator ξ α ′ β ′ γ ′ which is antisymmetric with respect to β ′ and γ ′ we obtain 0. Thus, the Maxwell curl equation reduces to:

ξ α ′ β ′ γ ′ Λ α β ′ Λ β γ ′ ∂ α H β = -iωεE α ′ (8.133)
let us multiply both sides by Λ γ α ′ and make summation on dummy index α ′ . We obtain: we then obtain a set of equations relating the complex amplitudes of the field components where the F α and the G α play a fully symmetric role:

ξ γαβ det Λ α β ′ ∂ α H β = -iωεΛ γ α ′ E α ′ = -iωεE γ = -
ξ αβ γ ∂ β F γ = k √ gg αβ G α ξ αβ γ ∂ β G γ = k √ gg αβ F α (8.138)
where k = ω √ µε

Appendix 8.C: Summary of tensorial approach

In curvilinear coodinates systems, the Maxwell's equations are based on the tensorial formalism deduced from relativity. If we consider only materials which are stationary with respect to the coordinate system, then the four-dimensional formalism developed by Post can be simplified. Maxwell equations are written :

ξ αβ γ ∂ β E γ = -∂ t B α ξ αβ γ ∂ β H γ = ∂ t D α + J α ∂ α D α = ρ ∂ α B α = 0 (8.139)
Post's formalism preserves the affine nature of Maxwell' equations: their expression is independent of the coordinate system. The geometry only appears in the constitutive equations along with the material's properties

D α = ε αβ E β B α = µ αβ H β (8.140)
In a perfectly linear, isotropic media with permittivity ε and µ, these relation ships become:

ε αβ = ε √ gg αβ ε αβ = ε √ gg αβ (8.141)
where g αβ are the contravariant components of the metric tensor.

g αβ = g αβ Assuming a time dependence of the form exp(-iωt), in a source free region if we substitute the constitutive equations for the material8.141 into Maxwell equations in the covariant form 8.139, setting

F α = E α G α = -iZH α with Z = µ ε (8.145)
we then obtain a set of equations relating the complex amplitudes of the field components where the F α and the G α play a fully symmetric role:

ξ αβ γ ∂ β F γ = k √ gg αβ G α ξ αβ γ ∂ β G γ = k √ gg αβ F α (8.146)
where k = ω √ µε

Figure 8 . 1 :

 81 Figure 8.1: Geometry of the problem: Two media are separated by a cylindrical periodic surface, with period d 1 described by the equation z=a(x)

Figure 8 . 2 :

 82 Figure 8.2: Geometry of the diffraction problem. Sketch of the coefficients for scattering matrix

Figure 8 . 3 :

 83 Figure 8.3: Schematic representation of diffraction at two surfaces separated by a layer .

Figure 8 . 4 :

 84 Figure 8.4: Coordinate system in which coordinate lines are parallel to ∆ and to Ox axis

Figure 8 . 5 :

 85 Figure 8.5: Echelette grating in three different oblique coordinate systems Figure(8.5) shows three typical grating surfaces obtained with (8.103) and with φ = 90, 63.4349 and 40 respectively, the latter demonstrating the extreme overhanging forms possible for small φ without the double value problem implicit with Cartesian coordinates.

Figure 8 . 6 :

 86 Figure 8.6: Various parametric representations of a right angle triangular profile

Figure 8 . 7 : 1 Appendix 8 .

 8718 Figure 8.7: Comparison of speed of convergence for two parametric representation of a right angle triangular profile. Full line: x 1 0 = .5, η = 0, dashed line: x 1 0 = .4, η = .9. Other parameters are: θ = 25 • , λ = 1,h = d 1 = 1

  .118) vectors e α are tangent vectors along coordinate curve x α . The matrix formed by the coefficientΛ β ′ α

Figure 8 . 8 :

 88 Figure 8.8: Curvilinear Coordinates:covariant and contravariant components of a vector in a plane

  125)andg = det( g αβ ) = det(J ) 2 = det Λ α ′ β (8.126)The g αβ establish a connexion between the A α and the A βA α = e α . A β .e β = g αβ A β or A β = g β α A α (8.127) Appendix 8.B: Transformation of Maxwell's equations

  123, 231, 312 -1 for α ′ β ′ γ ′ = 321, 213, 132 0 otherwise (8.129) Then let us change the coordinates: ∂ β ′ = Λ α β ′ ∂ α and the components

- 1 g

 1 = (det)(g αβ )(8.142) In an arbitrary curvilinear coordinates systeme x α , if the surface separating two materials, denoted (1) and (2), coincides with a surface of coordinates x 3 = constant, for example, then the conditions of continuity are expressed quite simply tangential component continuity:a 1 (1) = a 1 (2) a 2 (1) = a 2 (2) (8.143)normal component continuity: a3 (1) = a 3 (2) (8.144)

  is the identity matrix. The combined S matrix of the top and bottom interfaces and of the layer is given by: S = S a j * S j, j+1 * S a j+1 = S a j * S j, j+1 * S a j+1 = S a j * S j, j+1 * S a j+1

	S 21 1	1 S 22 1	S 2 =	S 11 2 S 21 2	S 12 2 2 S 22	(8.86)
	The star product * is defined by:
	S = S 1 * S 2 -1 S 11 1 2 × S 21 2 S 22 2 S 22 1 I -S 11 1 + S 12 S 11 = S 11 1 S 12 = S 12 1 × I -S 11 1 2 -1 × S 12 S 21 = S 21 2 × I -S 22 1 S 11 2 1 -1 × S 21 S 22 = S 22 2 + S 21 2 × I -S 22 1 S 11 2 -1 × S 22 1 × S 12 2	(8.87) (8.88) (8.89) (8.90) (8.91)
	where I (8.92)
	and finally, it turns out that
	S 11 = S 11 a j + S 12 a j P ( j) + U 2 P ( j)-S 11 a j+1	(8.93)

S 12 = S 12 a j P ( j)+ U 2 S 12 a j+1 (8.94) S 21 = S 21 a j+1 U 1 P ( j)-S 21 a j (8.95) S 22 = S 22 a j+1 + S 21 a j+1 P ( j)-U 1 P ( j)+ S 22 a j (8.96)

  iωεg γβ E β(8.134)

	Finally we get: and	√ gg γβ E β √ gg γβ H β ξ γαβ ∂ α H β = -iωε ξ γαβ ∂ α E β = iωµ		(8.135) (8.136)
	setting	F α = E α	G α = iZH α	with Z =	µ ε	(8.137)


Thomas Durt

Estimate of the weight of empty space based on astronomical observations

Keywords: 95, 30, Sf -Gravitation astrophysics 42, 25, Hz -Interference, optical 45, 20, D--Newtonian mechanics

As a consequence of the equivalence principle and of the existence of a negative vacuum energy, we show how a weak but universal linear response of the vacuum to a local gravitational potential suffices to explain the main features of so-called Pioneer anomalies as well as the apparent departure from the third Kepler law which has been observed at the level of numerous galaxies.

Introduction. -It is nowadays widely accepted that it is impossible to explain the deviations from Kepler's third law exhibited by stars at the level of numerous galaxies without formulating extra-hypotheses. The two main hypotheses that are presently [1,2] competing are (1) the MOND theory and (2) the dark matter hypothesis. We propose a third approach (VIND for Vacuum Influenced Newtonian Dynamics) in which Newtonian dynamics is preserved while, in the absence of dark matter, it is the linear response of the vacuum energy which is supposedly responsible for apparent deviations to Kepler's third law. The paper is structured as follows. In the first section, we introduce the MOND, dark matter and (3) VIND models. In the second section we investigate the modifications to standard Newtonian dynamics deriving from our VIND model at small distances compared to the typical scale of the model (L). This allows us to estimate the value of L, on the basis of so-called Pioneer anomalies. We then investigate the modifications at large distances and obtain a good agreement between our predictions and MOND phenomenology. The last section is devoted to open questions, discussion and conclusion.

MOND, dark matter and VIND models. - (1) In the MOND theory [3] ("MOdified Newton Dynamics") one assumes that Newton's law must be modified in a ad hoc fashion whenever we consider very small accelerations (relatively yo the inertial frame attached to "fixed" stars). It is motivated by the observed violations of Kepler's third law, which establishes that, as is the case in the solar system, celestial objects follow trajectories for which "the square of the period T varies like the cube of the mean distance R to the sun". In the case of circular orbits, this reflects Hooke's law: 2π T 2 R = a N ewton is proportional to 1 R 2 , up to constants. In other words, the Newtonian force m • a undergone by a test-particle of mass m at distance R of a ponctual source varies like 1/r 2 . On the contrary, MOND theory in its simplest form predicts that the (intensity of the) acceleration undergone by a particle test in a very weak gravitational field (for instance far away from other gravitational sources) obeys

|a MON D | = a N ewton • a 0 , (1)
where a 0 is a new universal constant introduced by Milgrom. It is easy to check that, far away from a ponctual source of mass M , MOND theory predicts that

a = a MON D = -G N M • a 0 /r 2 = -( G N M • a 0 )/r.
(2) Therefore, the acceleration undergone by a test-particle at distance R of a ponctual source varies like 1/R so that ( 2π T ) 2 R is attractive and proportional to 1 R . In other words, the velocity R/T is constant (does not depend on R) for any circular orbit, which is a fundamental property of MOND theories, that has been sometimes infirmed in the practice but has nevertheles been confirmed by numerous observations [4,5] (this is sometimes called the flat velocity distribution).

(2) in the dark matter approach it is assumed that the departures from Kepler's third law are due to the presence of a huge concentration of dark matter, invisible to our telescopes and accelerators but responsible for the appearance of large accelerations which cannot be explained otherwise, if we stick to the Newtonian paradigm.

Both approaches present difficulties. The MOND approach can be considered as a clever and economic phenomenological description of what happens in a large number of galaxies, but it suffers from a clear theoretical background and motivation. For instance it is very difficult to conciliate MOND theories and general relativity, the latter predicting, roughly, that at slow velocities, Newtonian gravitation prevails. It also suffers from clear experimental confirmations at the scale of clusters of galaxies and at the scale of the solar system.

The introduction of dark matter can also be considered to be an ad hoc hypothesis being given that, so far, all attempts to find a trace left by this hypothetical dark matter outside from the specific context of galactic rotations curves failed. At CERN for instance, none exotic particle described by dark matter theorists has been observed so far.

(3) The goal of the present paper is to show that the MOND predictions (or at least an "average" version of them) can be infered from an assumption that has the advantage to respect Newton dynamics. The price to pay is that we must assume that, everywhere in the cosmos, empty space exerts a weak back-reaction onto the gravitational potential. Otherwise, our results are derived in the framework of Newton's theory of gravitation. In particular we do not need to modify Newton's dynamical law (F = ma) at small accelerations, and the gravitational potential still evolves according to Poisson equation

∆φ = 4πG N ρ, (3)
where φ is the gravitational potential, G N the Newton gravitational constant and ρ the density of mass-energy. It has been recognised since several years that in the Newtonian limit the density of mass-energy ought to contain a contribution due to the so-called cosmological constant, which leads to the modified Poisson equation [START_REF] Straumann | Séminaire Poincaré[END_REF] ∆φ = 4πG N (ρ matter -2ρ λ ),

where ρ λ = λ/8πG N and λ is Einstein's cosmological constant. It is not completely clear whether this contribution is due to quantum vacuum fluctuations or not but here we shall take for granted that empty space possesses a negative mass density, constant in space. The novelty of our approach is that we also assume that, as a consequence of the mass-energy equivalence principle, the matter and the vacuum do not only contribute to the density of mass-energy as isolated source terms, but that they also exert a back reaction proportional to the gravitational field:

∆φ = 4πG N (ρ matter -2ρ λ )(1 + φ/c 2 ) = 4πG N (ρ matter -2ρ λ ) + (4πG N c -2 ρ matter -k 2 )φ, , (5)
where

k 2 = 8πG N c -2 ρ λ = λc -2 .
We are not interested here in computing effects due to the linear in φ back reaction of matter (4πG N c -2 ρ matter φ). They can be evaluated independently and in general they only result in a small correction to usual Newtonian dynamics. For instance, at the level of the solar systems where the gravitational potential is everywhere small relatively to c 2 , they does not modify much the Newtonian, zero-order dynamics. Nevertheless, the back reaction due to the vacuum energy, despite of the fact that it is weak, is present everywhere in space and its cumulative effects may not be neglected. Our goal is to show how these effects provide a plausible explanation of the observed deviations to Kepler's third law, and of so-called Pioneer anomalies as well... In summary, Vacuum Influenced Newton Dynamics derives from the assumption that Newton dynamics is not modified but "only" indirectly influenced by the appearance of a new source term k 2 φ that we ultimately attribute to vacuum1 .

VIND in the near and far field regimes . -Formally, the equation ( 5), in vacuum, can be solved once we know the associated Green's function, which is a solution of Green's equation

∆G(x, x ′ ) + k 2 G(x, x ′ ) = δ(x -x ′ ), (6)
where δ(xx ′ ) represents the three-dimensional delta function. As is well-known, even if we impose that the Green function is isotropic (that is to say, it depends only on the distance between x and x ′ ), it is not unique in the sense that all Green functions differ by a solution of the homogeneous Helmholtz equation. We are interested in a real and isotropic Green function, which suggests the choice

G(x -x ′ ) = -cos(k(|x -x ′ |) 4π|x -x ′ | . (7)
We are free to add to the Green function defined by equation ( 8) a real multiple of the harmonic function

sin(k(|x-x ′ |) 4π|x-x ′ | .
However, in realistic situations, the weight of such a solution remains very small, and the essential conclusions of all our analysis remain unchanged.

In particular, convoluting the omnipresent vacuum energy -2ρ λ everywhere in space with the Green function yields a finite constant additive correction to the potential. The contribution in a location x essentially comes from a region situated at a finite number of wavelenghts (denoted L, L = 2π/k) around x. It is not necessary to make particular assumptions about the structure of the universe (e.g. is it finite, open or else) in order to renormalise properly the gravitational potential in the present case, excepted that L supposedly differs from 0. We shall thus consistently forget in what follows the presence of the source term -2ρ λ in equation ( 5) and focus on the changes brought by the linear vacuum response k2 φ.

Near field regime. At small distances (|xx ′ | << L), the Green function reduces at the zeroth order to the usual Green function of Newton's theory -1 4π|x-x ′ | , which means that our model is nearly equivalent to Newton's model for gravitation at small distances from the sources of matter. As it is granted that, in the solar system, Newton's model works very well 2 , we must impose that L is quite larger than the size of the solar system. It is however interesting to take account of slight but systematic deviations of Newton's theory which have been observed with slow objects (the Pioneer 10 and 11 probes), and constitute what is called the "Pioneer anomaly".

In short, it has been repeatedly observed at the level of these probes (twice) that a constant centripetal acceleration seemingly prevails in our solar system. This acceleration is constant in time and space and directed towards the sun. Its magnitude is [8] of the order of 8 • 10 -10 m/s 2 .

Our model precisely predicts such an effect, as can be seen by developing the Green function (8) to the lowest order in k which reads

G(x -x ′ ) ≈ -1 4π|x -x ′ | + k 2 8π |x -x ′ |. (8)
Sufficiently far from the sun (where xx ′ is in good approximation constant in magnitude and direction for all source terms inside the sun (such that |x ′ | < R sun )), so that we expect that a constant radial and centripetal force of intensity One light year being equal to 9, 461 • 10 15 m, we find that L is more or less one fifth of a light year. This justifies, a posteriori, why the Pioneer anomaly is so weak while the usual Newtonian theory is so efficient at the scale of the solar system. These results are in agreement with the constraints derived so far, which are:

G N • M sun • k 2 /2
(1) L differs from 0.

(2) L is quite larger than the size of the solar system. Actually, although a more detailed treatment is out of the scope of the present paper, it is not too difficult to convince oneself that, in the approximation where we treat the sun as a sphere where ρ matter is taken to be constant and equal to < ρ matter > sun , it is relatively easy to estimate the perturbation resulting from (the modification of the usual Newton-Poisson equation by) new source terms including the vacuum density, as well as the back reaction of matter and vacuum resulting from the equivalence principle (as described at the level of equation ( 5)). This refined model leads thus to the same prediction as before, which is, it predicts as main perturbation at the level of the solar system the existence of a constant acceleration directed towards the sun, of magnitude

G N • M sun • k 2 /2.
Coherent interferences in the far field regime.

If we consider a distribution of matter of size smaller than or comparable to the size of the solar system (like a typical star for instance), the potentials generated by each element of the matter distribution will coherently interfere in the "far field regime" (|xx ′ | >> L). A striking novelty of our approach, compared to the "standard" Newtonian approach is that, in this case, the gravitational field generated by a mass M located at the position x ′ , at the position x, is no longer the Newton-Coulomb potential

-GN M |x-x ′ | , but it is the oscillating function -GN Mcos(k(|x-x ′ |) |x-x ′ |
. This leads to several counterintuitive properties, the first of them being that gravitational forces are not necessarily attractive. Indeed, the acceleration a V IN D generated by a mass located at position x ′ onto a test-particle located at position x obeys

a V IN D = -∂ ∂r ( -G N M cos(k(|x -x ′ |) 4π|x -x ′ | ) = -kG N M sin(k(|x -x ′ |) |x -x ′ | - G N M cos(k(|x -x ′ |) |x -x ′ | 2 (10)
Therefore, in the far field regime (at distances larger than one light year), the dominating term of the potential created by any celestial object of size comparable to or smaller than the size of the solar system is, generically, the first one in the last part of equation (10):

a V IN D ≈ -kG N M sin(kr) r (11)
The expression (11) obviously differs from the expression (2). Now, as we have shown in the previous section, L being smaller than one light year, it is significatively smaller than typical galactic sizes (even smaller than typical sizes of galactic cores), which implies that the potential generated by the stars of the galactic core will not always constructively interfere. As we shall now show, it is possible to find a good agreement between MOND and VIND models in this case (incoherent interferences in the far field regime), provided we interpret the MOND theory in a statistical sense.

Incoherent interferences in the far field regime. Another counterintuitive feature of our model is that gravitational potentials generated by different sources of matter do not always add to each other in a simple way but have the ability to interfere with each other. This property is obvious if we consider electro-magnetic potentials but is totally absent from the Newtonian paradigm. To illustrate the consequences of this property, let us assume that the core of a galaxy is densely populated so that the law of large numbers is valid. The acceleration generated by the core of the galaxy onto a test particle (or another star) located in the suburbian zone of the galaxy will then obey (where the symbol ...≈... must be understood here as...is equal "in average", or is "typically" equal to...):

|a V IN D | ≈ -C N ⊙ G N k < M ⊙ > r , ( 12
)
where C is a constant of the order of 1 while < M ⊙ > is the average mass of a star of the galactic core and N ⊙ its number of stars. It must be so because the phases which characterize the potentials generated by the different local clusters formed by the N ⊙ stars of the core will add incoherently, so that the average "phase factor" < N⊙ i=1 sin(kr i ) > obtained after summing over their incoherent contributions will be Gaussianly distributed around zero, with a variance σ equal to C • N ⊙ , with C a number of the order of unity 3 , as predicted by the law of large numbers.

Consequently, a V IN D , the (averaged over space) value of the acceleration towards the core of the galaxy, predicted in the VIND approach, will be statistically distributed around the value 0, with a spread of the order of

σ•kGN <M⊙> r = C √ N⊙GN k<M⊙> r
, as announced at the level of equation ( 12).

We find in this case that the ratio between < a V IN D > and a D obeys

|a V IN D | a MON D ≈ Ck G N < M ⊙ > a 0 = k G N < M ⊙ > 2a 0 . (13)
Injecting in equation ( 13) the typical values M sun ≈ 2 • 10 30 kg, a 0 ≈ 1, 2 • 10 -10 ms -2 and k ≈ 3 • 10 -15 m -1 , we find that

|a V IN D | a MON D = a P ioneer a 0 < M ⊙ > M sun . ≈ 2, 5 < M ⊙ > M sun .
(14) This demonstrates how, in a statistical sense, MOND phenomenology can be derived from the assumptions that we originally formulated about vacuum energy, which are the basis of our VIND model. That this ratio is not exactly equal to one is not amazing for two reasons:

(1) MOND's acceleration merely provides an order of magnitude. Indeed, in the framework of MOND's phenomenology it is well-known that one free parameter must 3 From now on we shall make the natural choice C = σ(sin(kr)) ≈

1/2 = 2π 0 sin 2 xdx/2π)
be adjusted, for each galaxy, in order to take account of the mass of the galaxy.

(2) the ratio aP ioneer a0 which appears in equation ( 14) is likely to increase due to various sources of decoherence and to decrease due to the presence of dark positive matterenergy. For instance, due to the spread of k one should expect that the effective MOND acceleration diminishes with the distance to the center of the galaxy, due to decoherence at large distances.

A fine study of the calibration factors which are introduced galaxy by galaxy provide a possible test aimed at discriminating VIND and MOND theories. In particular the scaling in terms of <M⊙> Msun is not the same in both approaches 4 . Other tests are discussed in the last section.

Open questions, discussion and conclusion . -Other possible experimental tests.

We mentioned three possible tests which confirm or could confirm that our ideas are well-founded (Pioneer anomalies, flat galactic velocity distributions and fine study of the calibration parameters and velocity curves), for which no new brutto experimental data must be collected [4,5,8]. Other tests could be realized in the futureas we discuss now.

Although both theories predict that the direction of the acceleration undergone by a test mass (for instance a star in the suburbs of a galaxy) is directed towards the centre of the galaxy, the VIND and MOND approaches differ by an essential feature: the intensity of this acceleration slowly varies over space in the MOND model while, in the VIND model, it undergoes a strong relative change (it can even change its sign) at a scale L (supposedly of the order of a fifth of light year). This implies another main difference between both theories: MOND theory is per se a deterministic theory while, due to the random distribution of constructive and destructive interference effects over space, VIND predictions are per se statistical. It is only in average that both theories lead to comparable predictions. This suggests several approaches aimed at discriminating between both explanatory schemes. The first of them consists of checking the sign of the centripetal acceleration undergone by stars located at the edge of dense galaxies, but this is not an easy task for experimentalists because we are dealing here with tiny velocities and accelerations. In the case of nearly circular orbits around the galactic core, the radial velocities are particularly small and vibrations along the radial axis predicted by the VIND model should be particularly slow. Now, even if the results confirm that the sign of the centripetal acceleration is nearly always positive, this does not necessarily mean that the VIND approach must be rejected. The reason is that the distribution of the trajectories of stars, in the VIND approach, is not necessarily homogeneously spread over space. Indeed, it is reasonable to suppose that the zones where gravitation is repulsive (roughly regions where < sin(kr) > is negative) are unstable because there the test particle undergoes a repulsive gravitational interaction which can not be counterbalanced by the centrifugal force. It is therefore likely that the regions where sin(kr) is negative are sparsely populated (which is confirmed by the fact that a flat velocity distribution has been observed: when the gravitational force is repulsive the distance to the galactic core and the tangential velocity are a priori uncorrelated). The latter prediction, in turn, is likely to be submitted to experimental observations because it would imply that inside a given galaxy populated regions of extent L would coexist with empty regions of comparable size. Actually, in order to derive the distribution of trajectories of stars of a given galaxy, in the VIND approach, one must ultimately describe the full history of these trajectories which is a very difficult task. Nevertheless, it is very unlikely that this history conspires to reproduce the deterministic predictions of MOND theory and therefore a measure of statistical fluctuations of the velocities of stars situated at the same distance of the core of a galaxy should suffice in principle to discriminate between MOND and VIND predictions.

Remark that the VIND theory also predicts the appearance of repulsive barrers around a star, at distances separated by integer multiples of L, which could maybe be revealed by a direct study of accretion disks of stars present in media where a lot of interstellar gas is present. This would be the firt experimental evidence of negative (repulsive) gravitational interaction.

Another situation of interest occurs when the law of large numbers is not valid, for instance when a small number of stars or a small number of galactic clusters are present. This is interesting because it is known that in this regime the MOND approach does not fit well with observations. Let us consider for instance that three stars are aligned, and separated by a distance larger than L, and let us treat the star situated in the middle as a test particle. In the MOND approach it will undergo an acceleration G N (M 1 /r 2 13 -M 2 /r 2 23 ) • a 0 . In the VIND approach however it will undergo an acceleration equal to

| -GN M1sin(kr12) r12 + GN M2sin(kr13) r13
|. Obviously these predictions have few in common.

Conclusion.

Before we end this last section, it is worth noting that our model indirectly provides an estimate of the cosmological constant λ = c 2 • k 2 ≈ (6π) 2 • 10 -16 ≈ 4 • 10 -14 s -2 , to compare with the lower value of 10 -35 s -2 , which has been derived on the grounds of rather recent astronomical observations [9] (aimed at measuring indirect manifestations of the existence of a non null cosmological constant predicted in the framework of general relativity), and with the larger value of 10 87 s -2 which derives from a cut-off at the Planck scale [START_REF] Straumann | Séminaire Poincaré[END_REF]. Despite of the fact that it is not clear today, to anybody, what is the meaning of the cosmological constant 5 and how and why it should be related to quantum vacuum fluctuations and quantum vacuum energy, we were pleased to elaborate on the idea that a negative and constant "etheric" distribution of mass is omnipresent in the cosmos. As we have shown, the recognition of this fact has implications at the level of astronomy which cannot be ignored, and it is not necessary to resort to sophisticated theories (like general relativity, exotic particles and so on) in order to figure out what some of these implications could be.

The most amazing lesson that we can draw from our analysis is that (in the case where our theoretical building would appear to resist to facts) it allows us to infer "the weight of empty space". Indeed, it predicts for the vacuum density of mass -2ρ λ = -k 2 c 2 /4πG N a value close to -0, 9•10 -3 kg/m 3 , more or less MINUS one gram per cubic meter...Once we opened the Pandora box, many questions spontaneously arise, e.g., "is k constant throughout time ?", "to which extent does its existence influence the cosmologic evolution (expansion, galactogenesis, distribution of the cosmic background and so on) ?", "can its value be infered from the standard particle model ?"...As we see, black clouds in the blue sky of astrophysics possibly generate rain and rainbows. * * * As an outsider to the field who essentially built his culture of the topic through the reading of magazines, wiki sites (e.g. the site of John Baez dedicated to vacuum's energy (http://math.ucr.edu/home/baez/vacuum.html), the wiki sites respectively dedicated to the cosmological constant and the MOND theory) and so on, the author wishes to apologize in advance in the case where the bibliography would be incomplete.

  must exist, which adds to the usual Newtonian solar attraction. Identifying this force with the Pioneer anomaly, we derive an estimate of k: k 2 = 2a P ioneer /(G N • M sun ). Injecting in this relation the typical values M sun ≈ 2 • 10 30 kg, a P ioneer ≈ 8 • 10 -10 ms -2 and G N ≈ 7 • 10 -11 m 3 s -2 kg, we find k = 2π/L ≈ 3 • 10 -15 m -1 . This provides an estimate for L: L ≈ 2 • 10 +15 m (9)

See ref.[START_REF] Blanchet | [END_REF] for a similar angle of approach.

Up to the so-called Pioneer anomalies[8] that we discuss soon, and up to relativistic corrections which, anyhow, remain out of the scope of the present work, merely aiming at describing the motion of slowly moving objects relatively to "fixed stars".

MOND's acceleration scales like (< M ⊙ >) 1/2 , while VIND's acceleration scales like (< M ⊙ >)1 .

Averaging the acceleration expresssed in equation (10) over a spatial period, we see that gravitation vanishes. Rather ironically,