Tryfon Antonakakis

Fadi Baïda

Abderrahmane Belkhir

Kirill Cherednichenko

Shane Cooper

Richard Craster

Guillaume Demesy

John Desanto

Gérard Granet

Editor Evgeny Popov
email: e.popov@fresnel.frwww.fresnel.fr

Boris Gralak

Sébastien Guenneau

Daniel Maystre

André Nicolet

Brian Stout

Gérard Tayeb

Fréderic Zolla

Benjamin Vial

GratinGs: theory and numeric applications

The basic idea of the differential methods consisits in projecting the electromagnetic field on a set of basic functions in order to reduce Maxwell partial differential equations into a set of ordinary differential equations. When working in a Cartesian coordinates, the natural basis consists of exponentials, using the periodicity of the optogeometrical parameters. Diffraction by a single aperture requires working in the basis of cylindrical Bessel functions [7.1], while diffraction by an arbitrary-shaped single object requires vector spherical functions [7.2] as a basis.

The first studies using the differential method [7.3] appeared in the late 1960s, initiated by the birth of the computers. These studies concerned the modeling of diffusion of particles in nuclear potential by using the separation of variables of the radial Schrödinger equation. The method was called "optical method" due to the similarity between the Schrödinger and the Helmholtz equations. The first applications to grating diffraction appear in 1969 [7.4], but accurate and converging results required combining the differential method with conformal mapping techniques [7.5]. The classical differential theory as known nowadays was formulated in [7.6, 7.7]. One can find a detailed review on the classical differential method in [7.8] It appeared that the classical differential theory suffered from severe numerical problems in transverse magnetic (TM) polarization, as well as for deep gratings. The first breakthrough was made in the first half of the 1990s, by introducing orthonormalization of the differential equations during their integration [7.9] and followed later by the so-called Rmatrix or S-matrix propagating algorithms [7.10]. The second breakthrough improved considerably the convergence in TM polarization for lamellar gratings, by introducing the correct factorization rules (see further on), at first by chance [7.11] and after that using theoretical arguments [7.12], closely followed by a generalization to arbitrary profiles [7.13]. A detailed review can be found in [7.14].

Maxwell equations in the truncated Fourier space

Let us consider a structure with two-dimensional periodicity along the x-and y-axis (Fig. 7.1) with periods equal to d x and d y . The modulated (grating) region extends in z from z min to z max . Inside that region, for a given value of the vertical coordinate z, the permittivity ε and permeability µ are periodic functions in x and y that can be projected on exponential Fourier basis: . We shall deal with a monochromatic (wavelength λ) plane wave incident on the structure with a wavevector:

inc 0 0 0 k ( , , ) = α β -γ  (7.2)
with components related to the incident polar angle θ (between the incident direction and the grating normal) and azimuthal angle ϕ (between the plane of incidence and the xOz-plane):

0 0 0 0 2 2 2 2 0 0 inc 0 0 0 k sin cos , k sin cos , k n , k 2 / α = θ ϕ β = θ ϕ γ = -α -β = π λ (7.3)
where n inc is the refractive index of the cladding. The existence and uniqueness of the solution of the diffraction problem is an interesting problem that is not discussed here. The reader can refer to several basic works (see for example [7.15, 7.16]. What is important to conclude is that the electromagnetic field is pseudo-periodic, so that similarly to eq.(7.1), the electric E  and magnetic H  field vectors can be represented in pseudo-Fourier series: In what follows, we use the notations: m 0 x n 0 y mK , nK α =α + β =β + . (7.5) From a numerical point of view, it is necessary to truncate the series in eqs.(7.1) and (7.4), introducing truncation parameters N x and N y , which limit the lower and the upper boundaries in the series.

Maxwell equations written in Fourier space take the form, assuming exp( i t) -ω time dependence with circular frequency ω: However, there are several problems in the numerical application of this rule:

i E (z) E (z) i B (z) dz d E (z) i E (z) i B (z) dz i E (z) i E (z) i B (z) d i H (z) H (z) i D (z) dz d H (z) i H (z) i D (z) dz i β - = ω -α = ω α -β = ω β - = -ω -α = -ω m, z,m,n m,n z,m,n y,m,n x,m,n x,m,n m,n z,m,n y,m,n d
First, numerical applications are simplified when using matrix notations. However, most of the standard routines use single-rank vectors and rectangular (2-ranks) matrices, while the vectors D and E in eq. (7.8) have two indexes, and the matrix ε depend on four indexes. In the case of classical grating with one-dimensional periodicity, this problem does not exist. Fortunately, for structures having 2D periodicity, a reduction to standard arrays is possible by introduction of a single index instead of the double for the vectors, by the following substitution:

x y y p (m N )(2N 1) (n N 1) = + + + + + (7.9) so that when m varies between -N x and +N x and n varies between -N y and +N y , p varies between 1 and P max = (2N x +1)(2N y +1). Using these notations, we can introduce standard arrays in the following manner: That can be written in matrix notations in the form

  D(z) (z) E(z)     = ε       , (7.12)
where double square brackets stand for the Toepliz matrix.

In addition, two diagonal matrices are useful: with p,p′ δ being the Kronecker's symbol. Second, due to the vectorial character of the fields, the matrix form in eq.( 7.12) has to be interpreted in a block form:

[ ] [ ]   [ ] [ ] x x y y z z D (z) E (z) D (z) (z) E (z) D (z) E (z)             = ε                     , isotropic media (7.14) [ ] [ ]         [ ] [ ]
D (z) (z) (z) (z) E (z) D (z) E (z) (z) (z) (z)   ε ε ε                 = ε ε ε                     ε ε ε                                    
, anisotropic media. (7.15) The third problem linked with the truncation of eq. ( 7.8) has limited the use of the differential methods (including RCW method) for more than 30 years, and has been solved for lamellar gratings in the late 90s [7.11, 7.12], and for arbitrary-profile gratings in the start of the 2000s [7.13]. The problem is due to the very slow convergence with respect to the number of Fourier components in the truncated sum of eq. (7.8), when the two functions in the product are discontinuous. As demonstrated by Li [7.12], four different cases can be distinguished with respect to eq.(7.12):

1. Both ε and E are continuous functions of x and y. 2. ε is discontinuous, but E is continuous. This is the case of the tangential component of E. 3. Both ε and E are discontinuous, but their product D is continuous, as it happens for the normal component of D. 4. All three functions are discontinuous.

In the first and second case, Laurent's rule assures relatively rapid convergence. In the third case, more rapidly converging scheme can be obtained through the following considerations for isotropic media.

If D = εE is continuous, then it is possible to factorize the product between D (continuous) and 1/ε (discontinuous) using the Laurent's rule (called by Li direct rule):

1 E(z) D(z) (z)     =     ε             , ( 7.16)
wherefrom the so called inverse rule is formulated:

1 1 D(z) E(z) (z) -     =     ε             , ( 7.17)
which can be applied if the matrix

1 (z) ε          
is not singular, a requirement that can create numerical problem for highly conducting gratings having small imaginary part of ε. (after [7.13]). When infinite series are considered, eq.(7.17) is identical with eq.(7.12). However, as shown in Fig. 7.2, the correct use of the direct or the inverse rules improves drastically the convergence of the differential methods with respect to the truncation parameter. Similarly to the abbreviation FFT, standing for Fast Fourier transformation, we have introduced the term Fast Fourier factorization (FFF) to name the correct use of the direct and the inverse rules, when applied numerically in the truncated Fourier space.

In the fourth case, neither the direct, nor the inverse rule result in acceptable convergence, so that this case must be avoided. Fortunately, this can be done by considering separately the electromagnetic field components, tangential or normal to the grating profile and taking into account that the electric field components tangential to the surface separating two different permittivities are continuous, in the same way as the normal components of the displacement D.

Differential theory for crossed gratings made of isotropic materials

In the isotropic case, the displacement vector D  can easily be separated into a continuous part

N N D E = ε  
, normal to the profile surface, and

T T D E = ε   that contains the continuous function T E  .
Let us define a unit vector N  , normal to the grating profile. Although it is well defined on the profile (except edges), it is necessary to generalize its definition all over the grating region, which cannot be done in a unique manner. Different choices are shown further on for specific gratings having 1D or 2D periodicity. Using this generalized vector, the relations between E  and D  can be decomposed into two terms, for each of which we can apply the direct or the inverse factorization rules, skipping the explicit writing of the z-dependence:

( ) ( ) N T D E E N N E E N N E   = ε + ε = ε ⋅ + ε - ⋅             . (7.18)
The first term is a product of type 2 and requires the direct rule. The second term is of type 3, demanding the inverse rule, so that:

    ( ) ( ) 1 T N 1 1 D E E 1 E N N.E N N.E . - -       = ε +       ε     = ε - +     ε                           (7.19)
Introducing a square matrix representing a tensor product denoted ( )

NN  
with elements given by N i N j , we obtain:

    1 1 D E NN E Q E - ε           = ε + -ε =             ε                       , ( 7.20)
where the matrix Q ε has the form:

    1 1 Q NN - ε   = ε + -ε     ε                   . (7.21)
In a similar manner for magnetic materials, we can find the link between magnetic field and induction in the truncated Fourier space:

    µ 1 µ B Q H , 1 with Q Q µ µ NN µ - ε     =         ==+ -                         (7.22)
Eq. (7.20) allows eliminating E z in the system (7.6):

[ ] [ ] [ ] ( ) [ ] [ ] 1 z ,zz z ,zx x ,zy y y x 1 ,zz ,zx x ,zy y E Q D Q E Q E H H Q Q E Q E - ε ε ε - ε ε ε   = - -       α -β       = - + +     ω   (7.23)
where the matrices α and β are defined in eq. (7.13).

Repeating the procedure for H z :

[ ] [ ] [ ] y x 1 z µ,zz µ,zx x µ,zy y E E H Q Q H Q H -     α -β       = - -     ω   , ( 7.24)
it is also eliminated from eqs. (7.6).

For non-magnetic media, the last expression is further simplified:

[ ] [ ] y x z 0 E E H µ   α -β   = ω . (7.25)
Thus the system (7.6) is replaced by a system of ordinary differential equations: .26) This equation can be expressed in a compressed form:

[ ] [ ] [ ] [ ] x x y y x x y y E E E E d iM dz H H H H                     =                         . ( 7
d F(z) iM(z)F(z) dz = (7.27)
Here the matrix M has 4x4 blocks:

1 1 ,zz ,zx ,yz µ,zz M Q Q Q Q - - ε ε µ = -α - β 1 1 ,zz ,zy ,yz µ,zz M Q Q Q Q - - ε ε µ = -α + α 1 1 ,yz ,zz ,zx ,zz ,yx M Q Q Q Q Q - - µ µ µ ε µ α = -ω + β + ω ω 1 1 µ,yz µ,zz µ,zy ,zz µ,yy M Q Q Q Q Q - - ε α = -ω - α + ω ω 1 1 ,zz ,zx ,xz µ,zz M Q Q Q Q - - ε ε µ = -β + β 1 1 ,zz ,zy ,xz µ,zz M Q Q Q Q - - ε ε µ = -β - α 1 1 µ,xz µ,zz µ,zx ,zz µ,xx M Q Q Q Q Q - - ε β = ω + β -ω ω 1 1 µ,xz µ,zz µ,zy ,zz µ,xy M Q Q Q Q Q - - ε β = ω - α -ω ω (7.28) 1 1 ,yz ,zz ,zx µ,zz ,yx M Q Q Q Q Q - - ε ε ε ε α = ω - β -ω ω 1 1 ,yz ,zz ,zy µ,zz ,yy M Q Q Q Q Q - - ε ε ε ε α = ω + α -ω ω 1 1 µ,zz µ,zx ,yz ,zz M Q Q Q Q - - ε ε = -α - β 1 1 µ,zz µ,zy ,yz ,zz M Q Q Q Q - - ε ε = -α + α 1 1 ,xz ,zz ,zx µ,zz ,xx M Q Q Q Q Q - - ε ε ε ε β = -ω - β + ω ω 1 1 ,xz ,zz ,zy µ,zz ,xy M Q Q Q Q Q - - ε ε ε ε β = -ω + α + ω ω 1 1 µ,zz µ,zx ,xz ,zz M Q Q Q Q - - ε ε = -β + β 1 1 µ,zz µ,zy ,xz ,zz M Q Q Q Q - - ε ε = -β - α .
This form looks like the form of the M-matrix obtained by Lifeng Li for crossed anisotropic (electrically and magnetically) gratings with profiles invariant with respect to z [7.17]. Whatever the form of the matrix M, eq.(7.26) represents a linear set of first-order ordinary differential equations. It can be solved numerically (with several problems, discussed further on), using well developed numerical schemes. In the case of vertical invariance of the optogeometrical parameters of the system inside the modulated region, the elements of the Mmatrix becomes constant in z, so that the solution of eq. ( 7.26) can be found through the eigenvectors and eigenvalues of M, a technique known under the name of Fourier modal method, or Rigorous coupled wave (RCW) method.

The solution of (7.26) gives a linear link between the field in the substrate and in the cladding max min

F(z ) T F(z ) = , (7.29)
where T is called transmission matrix. The advantage of this presentation comes from the fact that the field components participating in the calculations are tangential to the interfaces between the substrate and the modulated region, and between the cladding and the modulated region, so that they are continuous across these interfaces (in the absence of surface charges).

Electromagnetic field in the homogeneous regions -plane wave expansion

In most case, the substrate and cladding are homogeneous isotropic media. The electromagnetic field there can be expressed as a sum of plane waves. In particular, if the x and y-dependencies are given as in eq.(7.7), the z-dependence is explicitly known, for example for the electric field it takes the form: p p p p p E (z) A exp(i z) A exp( i z)

+ - = γ + -γ    (7.30)
of two waves propagating upwards (sign +) and downwards (sign -) along the z-axis, with p given in eq.(7.9). Each diffraction order with a given p propagates independently of the others, the coupling is effective inside the grating region. The z-propagation constant γ depends on the medium properties:

1 H E E 1 H E E   α β β + γ   = - +   ±γ ωµ ωµ     α + γ α β   = +   ±γ ωµ ωµ   (7.32)
where the sign of γ determines the direction of propagation in along z-axis. With this link in mind, the column vector F in eq.(7.27) takes the form:

[ ] [ ] x y x y E E F A A H H + + --           ≡ = Ψ + Ψ             , ( 7.33)
where the column vectors

x y A A A ± ± ±         =           (7.34)
contains the amplitudes of E x and E y propagating in positive or negative direction of the zaxis, matrices ± Ψ are block-diagonal:

xx xy yx yy , ± ± ± ± ±       Ψ =   Ψ Ψ     Ψ Ψ     (7.
± ± ± ± = α β β + γ Ψ = Ψ = γ ωµ γ ωµ α + γ α β Ψ = ± Ψ = ± γ ωµ γ ωµ    (7.36)
found from eq.( 7.32) Let us consider the case of a single incident wave from the cladding. The grating generates different diffraction order that propagate upwards in the cladding and downwards in the substrate. We attribute number 1 to the substrate and number 3 to the cladding. The total number of unknown diffracted field amplitudes will be equal to 4P max , two sets of 1 x,p A -and 1 y,p A -transmitted in the substrate, and two sets of 3 x,p A + and 3 y,p A + . These unknown amplitudes are subjected to 4P max number of linear algebraic equation in (7.29).

In order to obtain the T-matrix, the numerical integration of eq.(7.26) is made by using the so-called shooting method, which consists of choosing 2P max linearly independent representatives of the transmitted field. These representatives must correctly reflect the link between the electric and magnetic field components, as given by eqs.(7.32). A typical example for the shooting vectors starting from the substrate is that matrix 1+ Ψ , which has 2P max linearly independent columns. Here again, the number 1 indicates the substrate.

Thus the F column vector at z = z min can be formally written as a linear combination of the unknown amplitudes A 1-:

1 1 min F(z ) A -- = Ψ  , (7.37)
Assuming that there is no incidence from the substrate side. Here the tilde indicates that the vector F is not yet the true solution of the diffraction problem. The result of the numerical integration from z min to z max will provide the values of F  at z = z max , which are also a linear combination

1 max F(z )A - 
of A 1-, due to the linearity of the problem. On the other side, the column vector F at the upper interface is equal to

3 3 3 3 A A + + - - ψ + ψ
, according eq.(7.33), thus a linear set of algebraic equations for the unknown amplitudes 1 A -and 3 A + is obtained, with the free part determined by the wave incident from the cladding side:

3 3 3 3 1 max A A F(z )A + + - - + ψ + ψ =  . (7.38)
Once this system is solved, all field components can be calculated.

Unfortunately, this simple procedure creates enormous numerical problems that can be explained by using two different arguments:

First, it is known (but no quite well) in the theory of systems of ordinary differential equations, that numerical integration could become instable after a specific integration length, due to the fact that the set of shooting vectors can lose its linear independence during the integration. In other words, if the initial choice covers a vector space of 2P max dimensions, this space could shrink during the numerical integration to reduce its dimensions, so that the final algebraic system (7.38) could become singular. A solution of the problem based on this understanding was proposed in 1990 by G. Tayeb by using intermediate orthonormalization procedures during the numerical integration.

The second argument is based on the fact that inside the modulated region, as well as in the homogeneous regions, electromagnetic field contains components that propagate both in the positive and in the negative z-direction. During the integration, they both are treated in the same manner. As far as the solution requires taking into account the evanescent orders in addition to the propagating ones, a part of the former grows exponentially in z-direction, while the other part decreases exponentially. Due to the limited length of computer words, the ones that decrease substantially will be lost with respect to the ones that grow rapidly, even if the former could bring physical information. During the 90s, several different algorithms were proposed for solving the problem, based on a different treatment of the diffraction orders propagating upwards and downwards [7.9, 7.10]. Among them, to so called S-matrix propagation algorithm [7.10c] is probably the easiest to implement. Moreover, it can be used with methods other than the differential one in, for example, treating a stack of layers by the integral method, or by methods based on a transformation of the coordinate system. Interested reader can find in Appendix 7.1 a brief description of the S-matrix algorithm.

Several simpler isotropic cases

In practice most applications use non-magnetic materials, for which the form of M-matrix is considerably simplified, taking into account that then Q µ is diagonal and equal to 0 µ . Furthermore, several specific cases are of great interest for application, and they lead to a further simplification of the M-matrix.

Classical grating with one-dimensional periodicity, example of a sinusoidal profile

Let us consider a classical grating with grooves parallel to the y-axis and surface profile given by the equation z = g(x). The vector normal to the surface is given by ( )

2 1 N g (x), 0, 1 , if g (x) exists, 1 g (x) N (1, 0, 0), if not ′ ′ = - ′ + =   . (7.39)
where the prime stands for a derivative with respect to x. In case of vertical walls N (1, 0, 0). =  Thus the easiest way to generalize the normal vector to the entire modulated region is just to make it equal to eq.(7.39) not only on the profile z = g(x), but everywhere inside the grating region for min[g(x)] z max[g(x)] ≤ ≤ . The advantage of this choice is that N  does not depend on z, and the Fourier transformation of the tensor NN   is done only once. If the derivative of the profile function does not exist, or if the function is a multivalued one (e.g., circular or elliptical rods), but the interface can be expressed as a two-variable function: g(x, y) 0 = , (7.40) the normal vector is easily defined as the gradient of the profile function

[ ] [ ] N grad g(x, z) grad g(x, z) =  .
The Q ε matrix takes the form:

              1 1 2 2 z x x z 1 1 2 2 x z x z 1 1 N N 0 N N Q 0 0 1 1 N N 0 N N - - ε - -       ε + -ε     ε ε       = ε         -ε ε +       ε ε                                                             (7.41)
where it is taken into account that 2

2 x z N N 1 + = .
The fact that the normal vector components participate in the form of products in couples is important, because it leads to the conclusion is that the choice of the sign of N  plays no role. Further simplification of the M-matrix comes if limited to non-conical diffraction with

β 0 = 0:   1 1 ,zz ,zx ,zz 0 0 0 1 1 ,xz ,zz ,zx ,xx ,xz ,zz Q Q 0 0 Q µ 0 0 µ 0 M 0 0 0 µ Q Q Q Q 0 0 Q Q - - ε ε ε - - ε ε ε ε ε ε α   -α - α + ω   ω   -ω   =   αα -ω ε   ω       -ω + ω - α     (7.42)
This shows that the system to integrate decouples into two subsystems, corresponding to the two fundamental polarizations, transversal with respect to the plane of incidence, transverse electric (TE):

[ ] [ ]   y 0 x 2 x y 0 d E i µ H dz d H i E dz µ   = -ω     α   = -ω ε       ω   (7.43)
and transverse magnetic (TM):

[ ] [ ] ( ) [ ] 1 1 x ,zz ,zx x ,zz 0 y 1 1 y ,xx ,xz ,zz ,zx x ,xz ,zz y d E i Q Q E i Q µ H dz d H i Q Q Q Q E iQ Q H dz - - ε ε ε - - ε ε ε ε ε ε α     = -α - α -ω     ω       = ω - - α     (7.44)

Fourier transformation of the permittivity

The set of ordinary differential equations to be integrated contains the Fourier transforms of ε, 1/ε, µ, 1/µ, 2

x N , and 2 z N . In general, Fast Fourier transform (FFT) techniques can be easily applied. As already discussed with respect to eq.( 7.39), the normal vector components must be transformed only once, if chosen to be independent on z. On the other hand, the permittivity and permeability depend on z and their Fourier components have to be calculated for each value of z during the numerical integration. Fortunately, in the 1D case, it is possible and recommended to use analytical formulae for the Fourier transforms of ε, 1/ε, µ, 1/µ, which give faster more accurate results. This can be done because for a given value of z, they are piecewise constant functions of y. Fig. 7.3 presents schematically a grating with a period d that separates two homogeneous media with permittivities ε 1 and ε 3 . For a given value z 0 of z, the Fourier transform of, for example, the permittivity intside the modulated region 0 z h ≤ ≤ is given by

2 1 x x 1 2 1 2 x x d x imK x imK x 3 1 m x x 2 1 x x x imk 2 1 3 3 m,0 e dx e dx d d x x sin mK 2 ( ) e m + - - + - ε ε ε = + -       = ε -ε + ε δ π ∫ ∫ (7.45)
so that the two integrals can be solved analytically, once 1 x and 2 x are determined from the inverse of g(x):

1 1,2 0 x g (z ) - = . (7.46)
If the inverse of g(x) has more than two solutions, the sum of integrals (7.45) will contain several more terms. The same equations can be used to obtain the Fourier transforms of the inverse of the permittivity.

[ ] 0 x 1,2 2z h z 1 sin(K x) x arcsin( 1) 2 h = + ⇒ = -. (7.47)

Fourier transformation of the normal vector

As already explained, the Fourier transformation of the normal vector requires its continuation all over the space. If the grating profile can be represented as a single-value function, we can use eq.( 7.39) for N  and calculate the Fourier components of the tensor T NN   by use of the Fast Fourier transform (FFT) technique once for all z-values. For a sinusoidal grating having a profile defined in eq.( 7.47), the normal vector takes the form:

( ) x 2 2 2 x h cos(K x), 0,1 1 d N g (x), 0, 1 1 g (x) h 1 cos (K x) d π   -     ′ = - = ′ + π   +      (7.48)
.

Classical isotropic trapezoidal or triangular grating

A trapezoidal grating is shown schematically in Fig. 7.4 with two flat regions L at the top and the bottom of the groove and two different, in general, groove angles ψ. The Fourier transform of the permittivity and its inverse are calculated using eq.( 7.45) with:

1 0 1 2 C 0 2 x z cotg x x z cotg = ψ =- ψ (7.49) with C 2
x d L = -. For the normal vector, the period can be divided in four regions A to D, as shown in the figure:

y x 1 x 2 z 1 z 2 x z N 0 N sin N sin in A, in C, N cos N cos N 1 in B and D, N 0 = = ψ = ψ     = - ψ = ψ   =   =  (7.50)
Their Fourier components do not depend on z and can be represented as a sum of several analytical terms, similar to eq.(7.45): ( )

1 1 1 x x 1 2 x x 1 1 2 h cotg h cotg L 2 imK x imK x 2 1 x m 0 h cotg d L d 2 imK x imK x 2 h cotg L d L sin 1 N e dx e dx d d sin 1 e dx e dx d d ψ ψ + - - ψ - - - ψ + - ψ = + ψ + + ∫ ∫ ∫ ∫ (7.51) ( ) 1 2 x x 1 1 h cotg d L 2 2 imK x imK x 2 1 2 y m 0 h cotg L cos cos N e dx e dx d d ψ - - - ψ + ψ ψ = + ∫ ∫ (7.52) ( ) 1 2 x x 1 1 h cotg d L imK x imK x 1 2 x y m 0 h cotg L sin 2 sin 2 N N e dx e dx. 2d 2d ψ - - - ψ + ψ ψ = - + ∫ ∫ (7.53)
A triangular-groove grating can be considered as a particular case of a trapezoidal profile with no flat regions, 1

2 L L 0 = = , C x d = .
Moreover, the profile given in Fig. 7.4 also includes the case with vertical facets, and some more exotic profiles with hanging back walls, Fig. 7.5.

A B C D L 1 L 2 ψ 1 ψ 2 x 2 x 1 d x C

Classical lamellar grating

Lamellar profile with vertical walls is most easy to treat, because the normal to the profile vector has only one non-zero component, N x = 1. The Q ε matrix takes the form:

    1 1 0 0 Q 0 0 0 0 - ε     ε   =   ε     ε           (7.54)           1 1 0 0 2 1 1 0 0 0 2 0 0 0 1 2 0 0 0 0 0 0 µ 0 0 µ M 0 0 µ µ 1 0 0 µ µ - - - - - α α   ε β -ε α + ω   ω ω   β β   ε -ω - ε α   ω ω   =   α α - β -ω ε   ω ω     β β   ω - α   ε ω ω              (7.55)
In non-conical diffraction, when β 0 = 0, the two fundamental polarizations are decoupled and can be solved independently of each other. The M-matrix is simplified to obtain an antidiagonal block form:

    1 0 0 2 0 1 0 0 0 µ 0 0 µ 0 M 0 0 0 µ 1 0 0 0 - - α   -ε α + ω   ω   -ω     α =   -ω ε   ω     ω   ε             (7.56)
thus the two sets of differential equations for each polarization become:

[ ]   [ ] 1 x 0 y 1 y x d E i µ H dz d 1 H i E dz - - α     = ω -ε α     ω     = ω   ε          (7.57) and [ ] [ ]   y 0 x 2 x y 0 d E i µ H dz d H i E dz µ   = -ω     α   = -ω ε       ω   (7.58)
Even in the case of conical diffraction, it is possible to define two other polarizations, for which the differential system decouples. These are the electric and magnetic polarizations that are transverse with respect to the x-axis. Let us denote the two polarization with superscripts (e), when E x = 0, and (h), when H x = 0. For (e) case, it is possible to express H y as a function of H x from eq.( 7.26) and the first line of the M-matrix in eq.( 7.55):

    1 1 1 (e) (e) y 0 0 x H µ H - - - α α       =-ω -ε α ε β       ω ω    (7.59)
which can be simplified into:

  ( ) 1 (e) 2 2 (e) y 0 0 x H µ H -     = -α ω ε -α β     (7.60)
The second line of the matrix M then results in:

      ( ) 2 2 1 1 1 (e) 2 2 2 (e) 0 0 y 0 0 x d E i µ µ H dz - - -   β β     = ε -ω + ε α ω ε -α       ω ω    (7.61)
This expression can be further simplified, and together with the third line of eq.( 7.55) (when E x = 0) gives a set of equation for (e) polarization:

  ( )   ( ) 1 (e) 2 2 2 (e) y 0 0 0 x (e) 2 2 (e) x 0 y 0 d E i µ µ H dz d i H µ E dz µ -       = ω β ω ε -α -             = α -ω ε     ω  (7.62)
Similar procedure for (h) case when (h) x H 0 = , result in another system of differential equations, decoupled from the (e) case:

  ( )   ( ) 1 1 1 (h) 2 2 (h) y 0 0 x 1 (h) 2 (h) x 0 y d 1 H i µ E dz d i E µ H dz - - - -       = ω + β α ε α -ω       ε         = ω -α ε α     ω           (7.63)
In non-conical mount, β 0 = 0 and eqs.(7.62) and (7.63) become equivalent to eqs.(7.58) and (7.57).

Both conical and nonconical cases of diffraction by lamellar gratings are solved by eigenvector technique, due to the fact that the coefficients of the differential equations are zindependent. Moreover, due to the separation of the two fundamental polarizations, it is possible to further reduce by half the size of the matrices, by dealing with second-order differential equations. For example, eq. ( 7.57) can be written in the form:

[ ] [ ] x 14 y y 41 x d E iM H dz d H iM E dz   =     =   (7.64) Thus [ ] [ ] 2 x 14 41 x 2 2 y 41 14 y 2 d E M M E dz d H M M H dz = -     = -     (7.65)
Les us denote with 2 p ρ the eigenvalues of the product 14 41 M M and with V the matrix with its eigenvectors arranged in columns. The solution of the first eq.( 7.65) can be written as:

[ ] [ ] 1 x x E (z) V (z)V E (0) - = Φ (7.66) with pp pp p (z) exp( i z) ′ ′ Φ = δ ± ρ (7.
67) which shows that the elementary solutions along z (called modes, wherefrom the names Fourier modal method or Rigorous coupled waves method) exist in pairs that can propagate upwards or downwards with the same propagation constants.

By integrating the second eq.( 7.65), we obtain that:

1 y y H (z) W (z)W H (0) -     = Φ     (7.68)
with W that can be written in different forms, because the eigenvectors are defined within an arbitrary factor. For example, if we take into account the second eq.(7.64),

1 41 W i M V - = ρ 
, where the diagonal matrix ρ has elements equal to ρ p . Another possibility, that is quite convenient in TM polarization described by eq.(7.57), is at first to calculate the eigenvectors of 41 14 M M , instead of 14 41 M M (their eigenvalues are the same). Then the link between V and W contains the inverse of 41 M , which is just equal to

1 1 ω ε         , so that i 1 W V = ± ρ ω ε         .

Crossed grating having vertical walls made of isotropic material

Most of the recent applications of the Fourier modal method are devoted to studies of light diffraction by structures with 2D periodicity and piecewise invariant in the third direction. This popularity has several reasons. First, extraordinary light transmission was found in the late 90s by Ebbesen [7.18] for such structures, namely metallic sheets with periodical hole arrays, and it attracted a lot of attention (see Chapter 1). Second, the technology of such structures has significantly advanced in the last 20 years. Third, the Fourier modal method is relatively simple to implement, and much faster than most of the other methods, because of the eigenvalue/vector technique of integration.

Detailed study of these structures will be described in a separate chapter. However, due to its importance, we are discussing different aspects of this theory, as it presents a particular case of the more general geometry, that is characterized by a constant value of the zcomponent of the normal vector on each cross-section having z = const. The prolongation of the normal vector within the entire grating cell is discussed in sections 7.6.2.2.

Differential theory for anisotropic media

If we consider anisotropic media that do not extend inside the grating structure, there is not necessary to reformulate the diffraction theory, only that in the general case it is not possible to separate the problem into two independent polarizations, and it is necessary to work with the complete 4Pmax vectors and matrices.

In the case of anisotropic medium that lies inside the grating, the equations linking the M-matrix with the Q ε and Q µ matrices remain the same, eq.(7.28) for isotropic and anisotropic media. The difference comes from the fact that the Q-matrices take more complex form, because the link between the normal and tangential components of the couples E and D and H and B is made through the tensors of permittivity and permeability, which are not scalars. Let us establish this link in detail for E and D. As far as the continuous and discontinuous field components must be factorized differently, we construct a column vector F ε , which contains the continuous field components E T and D N . There are two tangential components to the grating surface, and only a single normal:

( )

1 2 N x T 1 y 2 z T N E D E F E T E U E T E E E ε ε     ⋅ ε           = = ⋅ =           ⋅                 (7.69)
where the double bar indicates a second-rank tensor with 3 dimensions, and the matrix U ε has the form:

( ) ( ) ( ) x y z 1x 1y 1z 2x 2y 2z N N N U T T T T T T ε   ε ε ε     =          (7.70)
with N ε  representing a tensor product with contraction of indices, for example,

( ) x xx y yx z zx x N N N N ε = ε + ε + ε  , etc.
The vectors 1 2 N,T ,and T    are defined on the grating surface, but for their further Fourier transform, it is necessary to choose a suitable continuation. The necessary conditions are that (i) they are continuous on the surfaces where ε and µ are discontinuous, and (ii) they form an orthonormal triad.

Since ε never vanishes, the determinant of U ε represents a quadric non-null form, equal to: ( ) ( )

1 2 i ij j i, j x,y,z det U N T T N N ε ε = ξ ≡ = ε ⋅ × = ε ∑    (7.71) since 1 2 N T T = ×    .
Thus U ε has an inverse inv U ε in the form:

( ) ( ) ( ) ( ) ( ) ( ) x 2 1 x x inv y 2 1 y y z 2 1 z z N N T N T 1 U N N T N T N N T N T ε ε       - ε × ε ×             = - ε × ε ×       ξ         - ε × ε ×                     (7.72)
It is not evident to derive this form, but it can easily be verified by using the equivalence inv U U ε ε =  and the fact that

1 2 N T T = ×   
. For example, the product of the second line of U ε with the second column of inv U ε can be written in vectorial form:

( ) ( ) ( ) ( ) ( ) ( ) inv 1 2 2 1 2 1 yy U U T N T N T T N T T N N ε ε ε ε     ξ =-⋅ ε × =- ε × ⋅ =-ε ⋅ × = ε ⋅ =ξ                (7.73)
Going back to the vector F ε , it is continuous across the grating surface, whereas the Cartesian components of the electric vector are, in general discontinuous, as well as the components of U ε . Thus for their Fourier transform, it is necessary to apply the inverse factorization rule:

[ ] 1 inv F U E - ε ε   =          (7.74) At the other hand, inv inv E U F D U F ε ε ε ε = ⇒ = ε   (7.75)
with F ε being continuous, so that the Fourier transform of D  requires the direct factorization rule:

1 inv inv D U U E - ε ε     = ε                (7.76) i.e., 1 inv inv Q U U - ε ε ε = ε          (7.77)
For gratings having anisotropic magnetic properties, the corresponding Q µ matrix is obtained from eqs. (7.71), (7.72), and (7.77) by replacing inv U ε by inv µ U and ε by µ .

Lamellar gratings made of anisotropic material

Such gratings are analyzed in the chapter devoted to the Fourir modal method by using more direct approaches, but here we want show how the corresponding equations can be obtained from the general eqs. (7.72). To this aim it is sufficient to realize that

N N ε × = ε -ε ε × = -ε ε ε ⋅ =ε       (7.79)
and eq.( 7.72) takes the form:

xy xz xx xx xx inv 1 U 0 1 0 0 0 1 ε ε   ε - -   ε ε ε       =           (7.80) with a determinant equal to xx 1 ε           . Thus 1 1 1 xy xz xx xx xx xx xx 1 inv 1 1 1 U 0 0 0 0 - - - - ε   ε ε   ε ε ε ε ε       =                                                                       (7.81)
The second matrix that is required takes the form:

1 0 0 U ε         ε ε ε ε ε   ε = ε - ε - ε ε ε     ε ε ε ε ε   ε - ε -   ε ε ε   (7.82)
This form is valid even when the permittivity tensor is not symmetric, as happens in the modeling of magnetooptical effects. The Q ε matrix takes the form obtained in [7.19], using a completely different approach:

1 inv inv Q U U - ε ε ε = ε          (7.83)
1 1 1 1 1 - - - - - ε ε ε ε ε ε ε ε ε ε ε ε ε ε = ε - + ε - ε ε ε ε ε ε ε                                                                                                                
- - - ε ε + ε ε ε ε ε ε ε ε ε ε ε ε - + ε - + ε ε ε ε ε ε ε ε                                                                                                        1 xz xx xx 1 -                 ε     ε ε                      

Normal vector prolongation for 2D periodicity; Fourier transform

As observed, the proper use of the direct and the inverse factorization rules requires that the vector normal to the interfaces between different media is defined not only on these interfaces, but throughout the entire grating cell. In the case of classical gratings with onedimensional periodicity, the prolongation of the normal vector can be done quite easily, as shown in sec.7.4.1. For two-dimensional periodicity, the choice depends on the geometry, but also on its mathematical representation. Several different solutions have to be considered, without pretending to be exhaustive.

In general, the cross-section profile changes with z, so that the matrices Q ε , Q µ , and M have to be recalculated for each value of z. If the geometry is z-invariant, this must be done only once. Concerning the Fourier components of the normal vector, there are two different classes of grating profiles that has to be treated separately. The first class consists of surfaces that can be expressed all over the unit cell (containing a single period in x and z) as an analytical (at lease piecewisely) function S z g(x, y) = , where S indicates that the point lies on the interface. In this case, it is possible to have a unique extension of N  whatever the values of z. In addition, it is not necessary to calculate the cross-section of the surface(s) with a plane perpendicular to the z-axis for each value of z. This case also includes multilayered homomorphous structure with constant layer(s) thickness in the z-direction. We consider this class of cases in sec.7.6.1.

The second class of surfaces includes surfaces that cannot be expressed through singlevalued functions, as the example given on the right-hand side of Fig. 7.12. In that case, it is necessary for each fixed value of z to know the cross-section function of the grating surface with the plane z = const. Subsection 7.6.2 presents general analysis, some important specific cases are considered further in the following subsections.

General analytical surfaces

If the interface representing the structure can be expressed as a single-valued function, analytical over the entire unit cell (this is also valid if different analytical functions can be defined over different regions of the cell): S z g(x, y), = (7.84) then the components of the vector normal to the surface defined on the surface have the form:

2 g g , ,1 x y N(x, y) . 1 grad g(x, y)   ∂ ∂ - -   ∂ ∂   = +  (7.85)
It can immediately be extended to whatever the values of z inside the modulated region. Moreover, its values do not depend on z, as it was a case of the classical one-dimensional grating already discussed in sec.7.4.1.

The permittivity and its inverse can easily be obtained on a mesh (x,y) covering the grating cell for each z:

3 1 g(x, y) z (x, y) g(x, y) z (x, y) < ⇒ ε = ε ≥ ⇒ ε = ε (7.86)
where 1 and 3 are the indexes representing respectively the inferior and the superior regions separated by the grating surface (as it was done in Fig. 7.3). If the cross-section of the grating surface with the planes z = const are ellipses (or circles), and if z N does not depend on (x,y) at each z, it is possible to replace the numerical Fourier transform by an analytical formulae. One important particular case is the z-invariant grating with elliptical cross section; another case includes the gratings having a rotational symmetry, as shown in Fig. 7.12.

The same extension (7.85) for the normal vector is valid for a stack of layers having homogeneous thicknesses in the z-direction: S, j j z g(x, y) t = + (7.87)

The permittivity and its inverse inside the intermediate layers are simply given as: S, j 1 S, j j z (x, y) z z (x, y) (x, y) .

-< ≤ ⇒ ε = ε (7.88)

Irregular general surfaces

If the case does not fit into the preceding section, the interface is expressed through the more general function S u(x, y, z ) 0 = , and the vector N  has to be determined for each inclusion:

S S S u u u , x y z N(x, y, z ) grad u(x, y, z )   ∂ ∂ ∂   ∂ ∂ ∂   =  (7.89)
However, these values are well defined on the grating surface (except on its edges), and have to be extended over the entire cell. When considerintg a cross-section of the profile with a plane at z = const., several different cases exist:

Single-valued radial cross-section

At first, we shall consider that the cross section function S S f (x , y ) 0

= defines a single curve, and that curve can be expressed in cylindrical coordinates as = ρ ϕ 5. Determining z N , together with N x and N y from eq.( 7.89).

Attributing these values of the components of S S

N(x , y )  to the pair (x,y). 7. Fast Fourier transform after the normal vector components are determined for all the pairs (x,y) on a mesh inside the grating cell. The procedure can be simplified for most of the typical diffracting objects, as shown further for objects with elliptical or circular cross-section.

If the grating profile varies with z, the calculations of the Fourier components of the permittivity and its inverse has to be made at each value of z, both for the analytical profiles, for which the normal vector prolongation can be chosen z-invariant, or for the irregular surfaces. For each (x,y) pair of the mesh used in the FFT method, it is possible to determine whether the point lies within or outside the cross-section part of Fig.

Objects with polygonal cross section

A typical example of such objects is presented in Fig. 7.1. Its surface consists of different plates, and for their treatment the condition N z = const. for fixed z is fulfilled, because N  is constant at each plate. Another possible surface consists of plane ribbons with curvature in zdirection, Fig. 7.7.

Fig.7.7. Surface made of plane ribbons

As shown in Fig. 7.8, the cross-section represents a polygon. On each of its sides the modulus of the in-plane component of the normal vector is known: assuming that N N =   i . The normal vector extension will be continuous everywhere, except on the sector border lines (bold dot-dashed lines) and thus the only points where both permittivity (and/or permeability) and N  are simultaneously discontinuous are at the polygonal corners, where anyway N  is never continuous. The procedure to follow requires that for each value of z the polygon corners (x i , y i ) and the z-components of N  for each segment are determined, as well as fixing some "central" point C C (x , y ) . Then the angular ranges of each segment with respect to that central point are calculated:

C C y y arctan x x - ϕ = - i i i (7.94)
For each pair (x,y), the azimuthal angle is given as

[ ] C C
arctan (y y ) / (x x ) ϕ = --, which value determines the number of the segment, say j, within the point lies. The unknown inplane part of the normal vector is perpendicular to the j-th segment:

( ) ( ) ( ) ( ) ( ) ( ) 2 z, j x, j j 1 j 2 2 j 1 j j 1 j 2 z, j y, j j 1 j 2 2 j 1 j j 1 j 1 N N y y y y x x 1 N N x x y y x x + + + + + + - = - - + - - = - - - + - (7.95)
The expression in the square root comes from the normalization of N  . The value of the permittivity depends on whether the point (x,y) lies inside or outside the polygon. The calculations of (x, y) ε and 1 ε are made simultaneously with the normal vector calculus. After the angular segment in which the point (x,y) of the mesh in grating cell is determined (say the j-th one, as in eq.(7.95)), we can find the length of S ρ between the central point and the polygon segment, shown in Fig. 7.8. For this sake we show in Fig. 7.9 the enlarged segment: Fig. 7.9. The j-segment of Fig. 7.

together with notations

The sine theorem gives the possibility to determine the angle ζ :

( ) ( ) ( ) ( ) 2 2 2 2 j 1 j j 1 j j 1 C j 1 C j 1 j
x x y y x x y y sin( ) sin( )

+ + + + + - + - - + - = ϕ -ϕ ζ (7.96)
wherefrom the radius S ρ is given as:

( ) ( ) 2 2 C j C j S j
x x y y sin( ) sin( ) Concerning the edges, in reality the surfaces never have such, as etching always ends by rounding the corners, as shown in Fig. 7.10. Let us consider that the rounding between the segments numbered i-1 and i is made preserving the values of N z , and that in the cross-plane z = const., the rounding can be considered as circular, having a center in C C (x , y )

- + - ρ = ζ π -ζ -ϕ + ϕ (7.
i i
and radius r i .The in-plane component of the normal vector follows the curvature radius and thus is given by equations, similar to eq.( 7.95):

( ) ( ) ( ) ( ) ( ) ( ) 2 z, x, C 2 2 C C 2 z, y, C 2 2 C C 1 N N x x y y x x 1 N N y y y y x x ∆ ∆ - = - - + - - = - - + - i i i i i i i i i i (7.99)
The prolongation is more complicated, if the consecutive values of N z at the two sides of the rounded corner differ significantly. In that case a linear interpolation of N z between ϕ i and ϕ + ∆ i i can be applied.

Mutlivalued cross-sections

If the cross-section cannot be represented as a radial function, nother possibility arises if it is a piecewise analytical function in x (or y), as shown in Fig. 7.11, where we can use two different functions of x. We assume again that z N is known, as it happens for z-independent profiles, for which it is simply null. In the upper part of the figure, for each value of x it is possible to calculate the normal vector on the profile:

( )

1 z 1 2 2 1 z f (x),1, N N 1 f (x) N ′ - = ′ + +  (7.100)
We can take this value to be the same for each y in the upper region 1 A , so that the numerical Fourier transform is made only once in 1 A and once in 2 A . The permittivity has to be calculated at each (x,y) mesh point:

S inside S outside y y , (x, y) y y , (x, y) .

< ε = ε ≥ ε = ε (7.101)

Objects with cylindrical symmetry

Many periodic systems consist of inclusion having rotational cylindrical symmetry, like spheres, vertical cylinders, or ellipsoids with axe of rotation parallel to the z-axis, but also smooth surfaces, as presented in Fig. 7.12. These structures are characterized by a circular cross-section of the surface with the horizontal planes at z = const., but also with an independence on x and y of the values of N z on each horizontal plane. In addition, due to the circular cross-sections, the angular component N ϕ = 0 everywhere. Once z is fixed, the variation of the interface in the vertical y x direction fixes the value of N z , for example through eq.( 7.89), wherefrom the radial normal vector component

2 z N 1 N ρ = - .
For each pair of (x,y) then: Fig. 7.12. Several profiles with cylindrical symmetry.

C x 2 2 C C C y 2 2 C C (x x )N N (x, y) (x x ) (y y ) (y y )N N (x, y) (x x ) (y y ) ρ ρ - = - + - - = - + - (7.102)
In addition, for profiles invariant in z-direction, N z = 0 and N 1 ρ = . The permittivity is given as a piecewise constant function:

2 2 2 inside C C 2 2 2 outside C C (x, y) , if (x x ) (y y ) R (z), (x, y) , if (x x ) (y y ) R (z) ε = ε - + - < ε = ε - + - ≥ (7.103)
where R(z) is the radius of the profile surface for a given z.

Having obtained the values of the normal vector components and permittivity for each x,y enables us to calculate their Fourier transforms, either by Fast Fourier transform (FFT), or analytically.

Objects with elliptical cross-section

Similar simplification is possible for systems with elliptical cross-sections that have N z = const. for z fixed. Such are the inclusions of vertical cylinders with elliptical cross-section, ellipsoids with one of the axes orientated in z-direction, but also all types of the structures shown schematically in Fig. 7.12 that have elliptical or circular cross-sections.

Let us assume that the ellipse axes are parallel to the x and y-axes.. The cross-section curve for z = const. is given by the equation:

2 2 2 S C S C x x y y R a b - -     + =         (7.104)
In order to obtain results similar to eq.( 7.102), we introduce an elliptical coordinates, defined as:

- = - -     +         - = - -     +         (7.106) with 2 z N 1 N ρ = -
, which remains constant for each fixed z. Concerning the permittivity, it is determined in the same way as in eq.( 7.103) for circular profile:

2 2 2 C C inside 2 2 2 C C outside x x y y (x, y) , if R (z), a b x x y y (x, y) , if R (z) a b - -     ε = ε + <         - -     ε = ε + ≥         (7.107)
with R(z) given by eq.( 7.104).

Remark on the prolongation of the normal vector

Special attention has recently been paid to the numerical implementation of the differential method for gratings having 2D periodicity formed by vertical holes or bumps that are invariant in z, and that have arbitrary cross-section in the xOy plane. A detailed study in the case of z-invariant geometry that applies for an eigenvalue/eigenvector technique of integration (FM or RCW method) is given in ref. [7.20], followed by several other works [7.21, 7.22]. It is necessary to note that the technique of prolongation of the normal vector as discussed in [7.20] can be applied also for z-dependent profiles with similar cross section; the difference is the renormalization factor 2 z

1 N for each z.

The authors compare several different formulations of the Fourier Modal method applied to structures with rectangular, circular, or elliptical cross-sections. These formulations include the classical formulation of Moharam and Gaylor that uses only the direct factorization rule, the formulation given by Lifeng Li [7.23] that introduces two different Fourier transforms of the permittivity ε, namely  ε        and  ε        , which are calculated by applying at first the inverse rule along one of the coordinates, and then the direct rule along the other one. This second formulation was made for rectangular and parallelogram crosssections. For circular or elliptical (or other smooth) forms, it introduces a stepwise treatment of the profile, which appears more slowly convergent than the special techniques developed after. The third approach to the problem requires a prolongation of the normal vector (NV) to the profile within the entire grating cell. As already stressed, there are several possibilities to make this. A typical example is the radial prolongation, Fig. 7.13a, which has been discussed in Sec.7.6.2.1 and 7.6.4 and it includes discontinuities of the normal vector on the cell boundaries, where the permittivity is continuous. Another approach proposed in [7.20] is the electrostatic one, which insures the continuity all over the cell and on its boundaries, except for on single points inside, Fig. 7.13b. Fig. 7.14 shows the convergence rates for the transmitted zeroth order of a grating consisting of dielectric cylindrical inclusions with a circular cross-section with refractive index n = 1.5, in normal incidence from the substrate. The grating period is 2λ, the width of the grooves is λ, and the grating depth is λ/(2n-1). The graph presents the diffraction efficiency in transmission as a function of the truncation order N using the three considered formulations: Moharam's original formulation, Li's formulation, and the formulation using the normal vector (NV) field. As usual the Fourier series run from -N to N, which yields 2N+1 Fourier coefficients for each of the two directions of periodic continuation, or 2 (2N 1) + coefficients in total. As can be expected, both the original approach and the formulation by Li have worse convergence than when correctly taking into account the factorization rules for the tangential electric field and normal displacement components to the profile, where the permittivity is discontinuous [7.20]. It is necessary to stress that the difference in the convergence rates is even more pronounced for metallic inclusions, having much larger optical contrast.

Fig.7.14. Convergence rates with respect to the truncation of the Fourier series for four different approaches used to model the diffraction by a cylindrical inclusion with circular cross-section.

Recently, Weiss et al. [7.24] proposed another alternative approach to treating smooth inclusions, by changing the coordinate system, so that its planes are parallel to the profile and to the grating sell walls (see Fig. 7.15). If the transformed system is orthonormal, its coordinate lines are automatically tangential or perpendicular to the physical walls. If not, the Maxwell equations have to be rewritten in covariant vector form using the covariant and contravariant vector components. Coordinate lines and surfaces according to (a) [7.20] and (b) [7.24] This approach is somehow equivalent to the normal vector prolongation, due to two main reasons:

(i)

The NV approach defines in an unambiguous manner the normal vectors on the profile, giving a liberty to continue them all over the cell. The coordinate transformation is also defined on the profile and the outside boundaries, but can be chosen in different ways around the grating cell.

(ii)

The change of the coordinate system introduces in the Maxwell equations the metric tensor  that multiplies the electric displacement and magnetic induction in the right-hand side of eq.(7.6), so that for the electric field we obtain the substitution:

D E E. = ε → ε     (7.108)
The normal vector approach acts in a similar manner by introducing the matrix Q , ε given in eq.7.21, which makes the following substitution in the Fourier space, eq.( 7.20

): D E Q E . ε       = ε →          (7.109)

Multiprofile surfaces

A grating with multiple bumps (or inclusions) inside the single cell could be treated by separation the cell into sub-cells, not necessarily rectangular, containing a single inclusion, as shown in Fig. 7. [START_REF] Wilcox | Scattering theory for the D'Alembert equation in exterior domains[END_REF], where a specific cross-section at z = const. is separated into three regions A, B, and C. As far as the Fourier components of the normal vector, of the permeability and the permittivity have to be calculated for each value of z (if they depend on z), the separation into subcells can vary with z. The case schematized in Fig. 7.17a can result from a surface covered by a thin layer of another substance, a layer that cannot be treated using eq.( 7 , where 1 N ( ) ϕ  is determined using one of the procedures discussed above for a single interface that is defined by the inner profile function.

If the point (x, y) lies in the outermost region, we take

2 N(x, y) N ( ) = ϕ  
, where 2 N ( ) ϕ  corresponds to the second interface. In-between, we have two possibilities. The first choice is to divide the region into two subregions as indicated in Fig. 7.17a with the dashed line. In each of them, N(x, y)  is taken to be equal to its values on the adjacent profile, so that it is continuous everywhere where the permittivity and/or permeability are discontinuous.

Structures with interpenetrating cross-section profiles

The second possibility is to introduce a linear interpolation inside region B, but it is necessary to know the distances 1 ρ and 2 ρ between the central point and the profiles along the ray with

[ ] C C arctan (y y ) / (x x ) ϕ = - - fixed. Then: 2 2 1 B 1 2 1 2 C C 1 2 N (x, y) N ( ) N ( ) N ( (x x ) ( ) , with y y ) ρ -ρ   = ϕ + ϕ - ϕ   ρ - + - -ρ ρ =     (7.110)
Another specific case that appears in the studies of magnetic resonators is presented in Fig. 7.17b. In can be treated in the same way is for the case in Fig. 7.17a, but it is necessary to introduce a separate region D indicated in the figure and containing the opening, for which N (0, 1, 0) =  , for example.

Integrating schemes

Numerical solution of a system of ordinary differential equations is a mature domain due to the enormous amount of physical and technical applications. Unfortunately, the grating problem represents one of the worst tasks for the theory of ordinary differential equations, because the system to be integrated is a stiff one. To better understand the problem, let us consider the case of a homogeneous layer that introduces no coupling between the diffraction orders. The solution of the diffraction problem contains waves propagating up-and downwards (in z-direction). These are plane waves, propagating or evanescent inside the layer. In lossless medium, their constant of propagation in z can be real, or imaginary, depending on the number of diffraction order under consideration: From the point of view of the theory of ordinary differential equations this means that the eigenvalues of the system differ significantly in magnitude, i.e., the differential system is stiff. The greater the difference, the more unstable the solution. On the other hand, the solution of the diffraction problem requires sufficient number of Fourier components of the profile function and electromagnetic field to be correctly represented by the truncated Fourier series, thus the necessity to work with large number of Fourier components, and thus the increasingly greater the stiffness of the differential system, i.e., more instable the solution with respect to the length and number of integration steps. The theory concludes that the so called explicit integration schemes are most instable for such problem, whatever their order, and implicit methods have to be used. The problem with the implicit methods is that they need one matrix inversion and several more matrix operations on each integration step, when compared with the explicit methods, so that the choice is not evident to ensure the most efficient integration scheme.

2 2 m 0 m (k n) γ = ± -α , ( 7
Let us recall the basic principle of the first-order explicit and implicit schemes. In a first-order approximation, the solution of the differential system: .114) between two consecutive points j z z = and j 1 z z + = can be searched in developing in series:

d F(z) M(z)F(z) dz = . ( 7
j 1 j j 1 j F(z ) F(z ) (z z )M(z)F(z) + + = + - . (7.115)
If M(z) and F(z) are evaluated in j z z = , this leads to the first-order explicit integration (Euler's) scheme:

j 1 j j F (z ) hM(z ) F(z )
+   = +    . (7.116)
where  is the unit matrix, and j

If M(z) and F(z) are evaluated in j 1

z z + = , we obtain the first-order implicit (inverted or backward Euler's) scheme:

1 j 1 j 1 j F(z ) hM(z ) F(z ) - + +   = -    . ( 7.117)
The theory says that this scheme is more stable, but it needs one matrix inversion on each step. A combination of the two must provide even better results, because it uses a half of the previous step:

1 j 1 j 1 j j 1 1 F(z ) hM(z ) hM(z ) F(z ) 2 2 - + +     = - +           . ( 7.118)
However, we need one additional matrix multiplication. In what follows we use these two single-point first-order methods under the names Expl 1 (single point explicit Euler integration) and Impl 1, eq.( 7.1 ) and compare the convergence with respect to the total number of integration points with several other more sophisticated integration schemes for two different metallic gratings in TM polarization, the most difficult combination when using the differential method.

The advantage of these formulations is that they all are single-step ones, and do not need a storage of the intermediate results on several integration steps. They can be easily programmed and don't need additional memory storage at each step. However, if we refer to one of the most relevant sources [7.25], we see that "this is the generic disease of stiff equations: We are required to follow the variation in the solution on the shortest length scale to maintain the stability of the integration, even though accuracy requirements allow a much larger stepsize." This means that a priori choice of the integration step without adaptive control and change in the step length cannot produce stable and relevant results. Unfortunately, it is quite difficult to use adaptive-step methods, because they require much longer computation times, as it is necessary to repeat the integration process several times when changing the integration step length. This is why we concentrate our attention to fixedstep algorithms.

Fixed-step multistep explicit methods have been used from decades in the differential method programming. The best results have been obtained when combined with an implicit correction by using a predictor-corrector scheme, as described further on. However, referring again to [7.25], "high order does not always mean high accuracy." It will be more useful, if larger integration step is obtained with high order or multistep methods, which is not obvious, as we observe on several numerical examples.

We have used three simple integration schemes, the single-point implicit or explicit scheme, as well as a 4-point predictor-corrector method (PCM 4). It contains two steps, the first one representing an Adams-Bashforth explicit 4-point scheme [7.26], described by the equation

+ + + + + + =   ′ ′ ′ ′ ′ + - + - +     . (7.119) with j j j F (z ) M(z )F(z ) ′ = .
(7.120) The corrector step is a 4-point Adams-Moulton integration:

j 4 j 3 j 4 j 3 j 2 j 1 j F(z ) F(z ) 251 646 264 106 19 h F (z ) F (z ) F (z ) F (z ) F (z )
, 720 720 720 360 720

+ + + + + + =   ′ ′ ′ ′ ′ + + - + -     (7.121)
which is an implicit scheme [7.27]. However, contrary to the other explicit schemes (BDF), it does not require inverting a matrix, it just makes one step back as a corrector. An extension of eq.( 7.1 ) to a multistep algorithm results in multistep implicit method, called also backward differentiation formulae (BDF) [7.28]. Typical 3-point and 5-point formulae take the form:

1 j 3 j 3 j 2 j 1 j 18 9 2 BDF3 : F(z ) hM(z ) F(z ) F(z ) F(z ) 11 11 11
- + + + +     = - - +        . ( 7
+ - + + + + +     = - - + - +        . (7.123)
It is evident that BDF3 (called further on Impl 3) requires a starting method for the first two points, and BDF5 (Impl 5), for the first 4 points. The same is valid for PCM in eqs.(7.1 ) -(7.1 ).

The second-order Runge-Kutta method is given in the form:

1 j j 2 j 1 2 j 1 j 1 j 2 k hM(z )F(z ) 1 k hM(z ) F(z ) k 2 F(z ) F(z ) k + + =   = +     = + (7.124)
which has an error proportional to 3 h . It is also called a midpoint point, because it requires the evaluation of the functions at the middle of the step, i.e., twice the number of the steps of the other tree methods discussed above.

The classical fourth-order Runge-Kutta method also uses midpoint values:

1 j j 2 j 1 2 j 1 3 j 1 2 j 2 4 j 1 2 j 3 j 1 j 1 2 3 4 k hM(z )F(z ) 1 k hM(z ) F(z ) k 2 1 k hM(z ) F(z ) k 2 k hM(z ) F(z ) k 1 1 1 1 F(z ) F(z ) k k k k 6 3 3 6 + + + + =   = +       = +       = +   = + + + + (7.125)
with an error of the order of 5 h . These two methods are higher-order explicit methods that do not need matrix inversions during the integration. However, as already stressed, higher order does not mean larger steps.

The highest possible order in h can be obtained theoretically, by using the eigentechnique:

( )

h 1 j 1 j F(z ) V e V F(z ) γ - + = (7.126)
where V is a matrix containing the eigenvectors of M and the exponential term in the round brackets represents a diagonal matrix constructed using the eigenvalues γ of M. In practice, this method does not increase the stability, because it remains a single-point explicit method. Moreover, it requires much longer computation time because of the requirement to solve an eigenvalue/eigenvector problem at each integration step.

The first example concerns a typical commercial sinusoidal aluminum grating that has very high efficiency in TM polarization. It supports a single diffracted order in -1 st order Littrow mount and has a modulation depth-to-period ratio of 40%. Fig. 7.20 presents a numeric test of the efficiency calculated for a different number of integration points using several integrations schemes. Due to the polarization and the grating material, it is necessary to separate the integration into several slices (5 in this case) in order to avoid numerical loss of precision, the results of the integration in two consecutive slices connected to each other by the use of the S-matrix algorithm. In Fig. 7.20b we have presented a part of the results, obtained with 20 slices, instead of 5. The comparison between the two cases show that 5 slices are sufficient, the weaker oscillation for 20 slices are due mainly to the fact that the horizontal scale is less dense, because the step in the total number of integration points is an integer times the number of slices. The truncation parameter N = 20, i.e., totally 41 Fourier harmonics of the field are used in the calculations.

As can be concluded, an absolute precision within 1% is rapidly obtained whatever the method used, with the total number of points of the order of 200. However, the predictor corrector method is less stable when the number of points is smaller than 300. Implicit methods are more stable, as expected, and result in an error smaller than 0.1% even for the number of points less than 200. It is interesting to observe than the first-order implicit method is more stable than the higher implicit methods, probably because it contains a middle-point evaluation of the field derivative, as seen in eq. (7.1 ). It requires a little bit longer computation time than the other two implicit methods, because of the additional matrix multiplication. The explicit method, which is the fastest one, shows slower convergence, as expected, whereas the performance of the higher-order RK methods competes with the implicit methods. Table 7.1 compares the computation times of the different methods for the two investigated cases (Fig. 7. and the following Fig. 7. ). For comparison, (null) indicates the time without any operation due to the integration, and that is necessary for the construction of the M-matrix and the use of the S-matrix propagation algorithm, as described in Appendix 7.A. The fastest method is the single-point explicit method, but as expected it is less precise given the same number of integration points, Fig. 7. . The implicit single-step middle-point method shows stability similar to the 4-th order Runge-Kutta method, but is slightly more rapid. The predictor-corrector method is less stable and requires longer computation times. It is necessary to stress out that in reality, the computation times are shorter than listed in the Table, because when the truncation N is smaller (usually 20 is sufficient), the number of slices for the S-matrix algorithm is smaller (due to the smaller number of evanescent orders taken into account); in addition, the total number of integration points used for constructing the Table are chosen to obtain 0.1% relative error, whereas in most of the cases just 1% is sufficient. The computation time grows linearly with the total number of integration points, as well as with the number of slices used in the S-matrix algorithm. The time dependence concerning the truncation N in the Fourier series grows as 3 N -3.5 N , because this parameter determines the size of the matrices.

When the total integration length is multiplied by 5, the number of integration points required is also multiplied by the same factor, as observed in Fig. 7. . A grating twice as deep as the period, acts almost like a flat mirror in TM polarization, with the efficiency in order -1 hardly exceeding 1%. Due to the large depth, the absorption is increased, so that the reflectivity in order 0 is equal to 56.78%. We compare the convergence in the weak -1 st order, so that even a small absolute error appears as a large relative error that can be easily observed in the figures. The number of Fourier harmonics (truncation parameter 2N+1) also has to be increased by a factor of 2.5 to 101 (N = 50). The number of slices in the S-matrix algorithm is increased seventh-fold to 35.

The first Fig. 7. compares several explicit integration schemes with the single-point implicit method. The main conclusion to be drawn is that the best scheme remains the implicit method, only the explicit Runge-Kutta fourth-order scheme seems to compete in convergence rate with respect to the total number of integration points, but somehow slower.

The comparison of several implicit methods confirms the general idea that multistep choice does not necessarily improve the stability (Fig. 7. ). When compared with Fig. 7. , the implicit methods are characterized by smaller oscillations when the number of steps is increased, but the most rapid convergence is obtained with the simplest procedure, single-step method (let us remark again that we use the middle-point calculations, as in eq.(7.1 )). Like all the other implicit methods, it requires a single matrix inversion on each integration step, but needs less memory storage, and avoids several matrix sums and multiplication by different constants, necessary for the multipoint methods. and apply than the theory of numerical methods of integration of ordinary differential equations.

The idea is sketched in Fig. 7.2 , where a sinusoidal surface-relief grating is approximated by a 5-stairs profile. While this approximation (with sufficient number of steps M, depending on the groove depth) works quite well in TE polarization, the TM case presents a convergence rate with respect the truncation number of Fourier components of the field much slower than the ordinary differential method (no staircase approximation), see Fig. 7.2 . Moreover, the greater the number of vertical slices M, the greater the truncation number required. A detailed analysis of this problem can be found in [7.14, 7.29], but the basic idea is quite simple. The staircase approximation substitutes the otherwise smooth sinusoidal profile by a profile that has sharp edges. The greater the number of stairs, the greater is the number of edges. It is well-known from general electromagnetism that edges introduce electric field singularities. While in TE polarization the only electric field components are tangential to the profile (in y-direction), thus have no discontinuities and singularities, this is not the case in At the edges of each step, a sharp maximum of the electric field is observed. These maxima are not a numerical artifact, they represent the physical effect of introducing edges to replace a smooth profile. These sharp variations of the field require larger number of Fourier components to be correctly represented. Moreover, the greater the number of slices (stairs), the greater the number of the maxima, thus the greater the truncation number required. Numerical experiment has shown that this phenomenon has nothing to do with the integration (eigenvalue/vector) technique, because the results of the convergence rate and field maps are the same for the staircase approximation when using the RCW technique or the differential method. On the contrary, if the true smooth profile is treated by the differential method by using a numerical integration of the ordinary differential system with the elements of the M-matrix depending on z, there is no such singularities of the electric field (Fig. 7.2 ), so that the convergence with respect to the number of Fourier harmonics is much faster, provide the correct factorization rules are used (Fig. 7.

).

Recently, some authors [7.30] have proposed to maintain the eigenvalue/vector technique, but to use the correctly determined Fourier presentation of the profile, i.e., the correct factorization rules, as presented in eqs. (7.39) -(7.44), instead of lamellar-profile factorization, eqs.(7.56) -(7.65), at each step. This is equivalent to using the formulation proposed by the differential method for a smooth profile, i.e. avoiding the field singularities at the edges, but to use the eigenvalue/vector technique of integration by assuming that the modified M-matrix, as given by eqs. (7.39) --(7.44), is z-invariant across each step height. We have already tried this in [7.29] and the conclusion was that using this approach, the number of steps (stairs) has to be relatively larger that by using some better adapted integration technique. And indeed, the eigenvalue/approach to a z-dependent system is equivalent to the rectangular rule with equidistant points of integration, one of the worst choices, as known from the theory of ordinary differential equations. In addition, due to eigenvalue/vector evaluation on each integration step, its computation times are several times longer than for the other methods (see Table 7.1 in the previous section), known from the theory of ordinary differential equations. This is why the authors of [7.30] need more than 2000 equidistant points of integration for a trapezoidal profile, for which the better adapted numerical integration scheme can suffice with 300 points. Unfortunately, the authors of [7.30] do not consider the differential method as a "reference method" in their work. Fig. 7.24. The same as in Fig. 7.23 but calculated using the differential method, after [7.29].

7.44

Gratings: Theory and Numeric Applications, 2012

Appendix 7.A: S-matrix propagation algorithm

Almost all electromagnetic theories work by providing the link between the electromagnetic field amplitude values established on two different interfaces. These values could be calculated in the real or the inverted space, or the projections of the field on some functional basis, etc. Whatever the theory, if the media are linear, the link can be expressed in a matrix form:

p p p 1 A T A - = . ( 7

.127)

Here, A stands for a column vector containing the field amplitudes in the given basis, the first interface has a number p-1, and the second on, p. p T is called transmission matrix between the interface (p-1) and p. Numerical problem arises due to the fact, that the "propagation" between different interfaces contains, in general, both growing and decreasing terms, due to both absorption losses or/and evanescent character of some field components. If a real field term propagates from p-1 to p (the green arrow in Fig. 7.A.1), it never grows (unless media with optical gains). Same is valid for the true propagation from interface p to p-1. However, eq. ( 7.1 ) is asymmetrical, i.e., it contains propagation only from interface p-1 to p, thus a naturally decreasing field that propagates in the opposite direction (from p to p-1), will be expressed in the T-matrix in the form of growing terms (the red arrow in Fig. 7.A.1). If the propagation length is sufficiently large, these artificial growing terms can overweight the other terms, mainly due to the finite numerical length of the computer word. One approach that overcomes this problem and that has become quite popular during the last 15 years is the so-called S-matrix propagation algorithm, S staying for 'scattering'. The basic idea is quite simple: As far as the problem of growing terms has been identified, let us try to do as Nature, by determining another link between the field amplitudes, by separating them into terms propagating (or decreasing) in direction (p-1  p) or in direction (p  p-1). Let us denote the first set with superscript +, and the second set by a superscript -. The S-matrix between the two interfaces provides the following link: The physical meaning is that amplitude p 1 A + -, which propagates from p-1 to p is defined on p-1 and is not growing in-between p-1 and p. In the same manner, the amplitude p A -that represents propagation from p to p-1 is defined on the interface p and is not growing in direction of interface p-1. To say in other words, the amplitude p 1

p p 1 p,p 1 p 1 p A A S A A + + - - - - -         =         . ( 7
A + -is incident on the interface p-1 from the previous interface p-2, the second amplitude p A -is incident on p from p+1, while the amplitudes on the left-hand side of eq.(7.1 ) are the amplitudes that are scattered in direction to the outside interfaces (p-2 and p+1), thus the name of the scattering matrix S. As observed in Fig. 7.A.2., the blocks p S --and p S ++ describe the physically correct transmission between p and p-1 or between p-1 and p, respectively, while the other two blocks, p S -+ and p S -+ describe the reflection on the interface p or p-1, respectively. This interpretation explains why there are no numerical problems due to the growing non-physical interactions when using the S-matrix. The advantage of this formalism is the absence of artificially growing terms in S. The inconvenience is that electromagnetic theories cannot give a direct expression of the matrix S. However, it is possible to express it by using the T-matrix elements, if it is possible to calculate them correctly. If the 'distance' between interface p-1 and p is quite large (with respect to the growing speed of the growing terms), there is loss of precision in determining the T-matrix. The problem can be solved by introducing additional artificial interfaces between p-1 and p in a such manner that to be able to correctly calculate the T-matrix in each subslice. Once the T-matrix calculated, the S-matrix can be obtained in a closed form. However, the total electromagnetic problem of diffraction (or scattering) requires the knowledge of the entire S-matrix of the system, because the physical problem to be solved needs to express the scattered fields as a function of the fields incident on the system (or generated inside, as is the case for electromagnetic antennas). There exists an iterative algorithm that enables us to establish the total S-matrix without calculating the elementary Smatrix between each consecutive pairs of interfaces, as stated in eq. (7.1 ). For that sake, we The initializing values of 0 S for p = 0 are just the elements of the unity matrix. As already said, it is necessary to be able to calculate the T-matrices for each intermediate medium between the interfaces. When advancing from the interface p to p+1, we obtain the T-matrix with subscript p+1: T -+ + can contain decreasing terms ('propagation from p to p+1), i.e., it could be numerically instable to invert them.

On the other hand, the 'next' S-matrix will link the amplitudes with index 0 to the amplitudes (p+1): At first, we express p 1 A - + from eq.(7.1 ) and substitute p A + from eq.(7.1 ):

( ) These relations exist in several possible forms, but this one is quite well adapted to the case without incident waves on interface 0, because in the iterative algorithm we need to calculate only the half of the blocks, namely the two given by eqs. (7.1 ) and (7.1 ).

p
The only matrix inversion in the iterative algorithm concerns the procedure to obtain the matrix  . The initial matrix 1 - contains the potentially large terms from p 1 T -- + , so that its inversion creates neither numerical problems to be inverted, nor growing terms to create numerical instabilities. ( ) The expression are quite similar in form to those obtained in Appendix A. Moreover, they allow avoiding the inversion of p 1 T +  . Finally, there exist a combination of expressions including partial T-matrices, treated separately in Appendix 7.A and 7.B. In some cases the link between the amplitudes on two consecutive interfaces or across a single interface that separates two different media can be expressed in the form: T . If this creates numerical problems (for thick layers), such direct approach is not applicable. In that case it is better advised to apply twice the S-matrix algorithm, at first in each direct form (Appendix 7.A), and then in the currently discussed inverted form. It is quite easy to understand the logic, by introducing a virtual set of amplitudes in eq.( 7.155):

0
p 1 p p 1 p 1 p A A A A + + + + - - +         =           T , ( 7

Fig. 7 . 1 .

 71 Fig.7.1. Schematical representation of a structure having two-dimensional periodicity in x and ydirections, consisting of truncated pyramids with height h.

  m,n z,m,n y,m,n x,m,n x,m,n m,n z,m,n y,m,n m,n y,m,n m,n x,m,n

6 )

 6 As can be observed, the third and the sixth equations are not differential equations, and they are used to eliminate the z-components of the fields, as shown further on. It has to be stressed out that the equations with different (m,n) numbers are coupled through the Fourier components of D E chapter we assume media with linear dielectric and magnetic properties and without spontaneous polarizations. The problem of Fourier transform of the product of two functions generally, solved theoretically by convolution of the Fourier transformers of the two functions, using the so-called Laurent's rule:

  ) D (z), E (z) E (z), etc. for H and B,

Fig. 7 . 2 .

 72 Fig.7.2.Convergence of the classical and the FFF version of the differential theory in the case of a dielectric sinusoidal grating with high contrast. Squares, old version of the differential theory for TM polarization; open triangles, new version, TM polarization; solid triangles, TE polarization(after[7.13]).

Fig. 7 . 3 .

 73 Fig.7.3.Piecewise constant representation of the permittivity for a one-dimensional grating

Fig. 7 . 4 .

 74 Fig.7.4. Trapezoidal profile with parameters. The normal vector direction is given in red arrows.

Fig. 7 . 5 .

 75 Fig.7.5. Two different profiles with slanted grooves

Fig. 7 . 6 .

 76 Fig.7.6. Single-curve cross-section of the grating surface at z = const.

7

 7

Fig. 7 . 8 .

 78 Fig.7.8. Object with a polygonal cross-section.

  its direction is perpendicular to the segment. If the i-th segment is located between the angles ϕ i and 1 + ϕ i , we can extend the definition of the normal vector all over the unit cell situated within the range ( ϕ i , 1 + ϕ i ), delimited by the bold dot-dashed lines in Fig.7.

Fig. 7 .

 7 Fig.7.10. Schematic presentation of corner rounding

Fig. 7 .

 7 Fig.7.11. Piecewise analytical c ross-section of the grating surface at z = const.

  these notations, the ellipse becomes a circle, for which the considerations of the previous subsection apply.

Fig. 7 .

 7 Fig.7.13. Two different prolongations of the normal vector for a circular inclusion. (a) Radial prolongation. (b) Electrostatic continuation of the normal vector for a circular cross-section inclusion inside a square grating cell (after[7.20]).

Fig. 7 .

 7 Fig.7.15. Coordinate lines and surfaces according to (a)[7.20] and (b)[7.24]

Fig. 7 .

 7 Fig.7.16. Cross-section at z = const. of a grating having different inclusions. The three different regions to be treated independently are separated by dashed lines.

1 N

 1 .87). The simplest possibility is to have different continuation of the normal vector inside each region. At first, the angle[] with coordinates (x, y) is calculated, and it is necessary to determine to which region the point belong. If it lies inside the innermost region C, the values of

  Fig.7.17. Structures with interpenetrating cross-section profiles

  of γ are bounded by 0 k n , but the imaginary parts are not bounded, as their asymptotic values for large |m| are given by:

Fig. 7 . 18 .

 718 Fig.7.[START_REF] Ebbesen | Extraordinary optical transmission through subwavelength hole arrays[END_REF]. Aluminum grating with period 0.5 µm and depth 0.2 µm used in -1 st order Littrow mount at 0.6328 µm wavelength in TM polarization. Convergence with respect to the total number of integration points, truncation to 41 Fourier harmonics and using 5 (a) and 20 (b) slices in the Smatrix algorithm. The acronyms for the methods are defined in the text.

Fig. 7 .

 7 Fig.7.[START_REF] Götz | Normal vector method for the RCWA with automated vector field generation[END_REF]. Schematical approximation of a sinusoidal grating profile, approximated by a 5-step staircase profile.(after[7.29]).

Fig. 7 .

 7 Fig.7.22. Convergence of the minus-first-order efficiency in TM polarization of the FMM (RCW)and the exact modal method (indicated on the figure) for a sinusoidal grating in a staircase presentation with M  = 20, as compared to the convergence of the differential method for a smooth sinusoidal profile (curve "diff."). Period d = 0.5 µm, groove depth a = 0.2 µm, aluminum refractive index n Al = 1.3 + i7.6, illuminated at 40° incidence with wavelength λ = 0.6328 µm, (after[7.29)].

  This can be observed in Fig.7.2 .

Fig. 7 .

 7 Fig.7.23. Spatial field distribution of |E x | 2 in the vicinity of several steps inside a groove of a 10step staircase profile, used to approximate the sinusoidal grating under study in TM polarization. The grating parameters are the same as in Fig.7.22, after [7.29].

Fig. 7 .

 7 Fig.7.A.1. Schematical representation of the action of the T-matrix between interfaces p-1 and p

Fig. 7 .

 7 Fig.7.A.2. Action of the S-matrix between interface p-1 and interface p.

  another intermediate S-matrix that corresponds to the scattering between some initial interface (numbered as 0) and the interface with number p

  7.1 )-(7.1 ) enable us to express the matrix p 1 S + as a function of p S and p 1 T + .

  case of the Fourier-modal (RCW) method across each interface, with the partial transmission matrices p 1+  Tcontaining the eigenvectors of the proper modes inside each media. The same expression is obtained in the coordinate transformation method when using eigenvalue technique of integration. Usually, in both approaches, one obtains the full transmission matrix by inverting p 1 +  T and multiplying the result by p 1 +

Table 7 .

 7 1. Computation times of the different methods described in the text for the two cases with groove depth to groove period equal to 40% and 200%

	Method	N = 20, slices 5	N = 50, slices 35
		int. points 400,	int.points 1500,
		modulation 40%	modulation 200%
	Expl. 1	1.41 s	72.3 s
	Impl.1	3.55 s	187.6 s
	Impl.3	2.67 s	157.8 s
	Impl.5	2.81 s	168.9 s
	PCM 4	2.14 s	120.0 s
	RK 2	2.70 s	144.9 s
	RK 4	4.70 s	252.5 s
	eigentechnique	7.54 s	360.0 s
	(null)	0.68 s	38.5 s

  It is obvious from the previous considerations that the growing terms are potentially present in the block p 1 T -- + ('antipropagation' from p to p+1), while the blocks p 1 T ++ + and p 1

			   	A A	p 1 p 1 + + -+	   	=	p 1 T +	   	A A	p p + -	   	.	(7.130)
	That will be expanded in the form:
	   	A A	p 1 p 1 + + -+	   	  =  	p 1 p 1 T T ++ + -+ +	p 1 p 1 T T +-+ --+	   	   	A A	p p + -	   	.	(7.131)

  7.28 A. Iserles, ed.: A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996 7.29. E. Popov, M. Nevière, B. Gralak, and G. Tayeb, " Staircase approximation validity for arbitrary-shaped gratings," J. Opt. Soc. Am. A 19, 33-42 (2002) 7.30. I. Gushchin and A. Tishchenko, "Fourier modal method for relief gratings with oblique boundary conditions," J. Opt. Soc. Am. A 27, 1575-1583 (2010)

											.156)
	 T	p 1 +	   	A A	p 1 p 1 + + -+	    =    	A A  	p p + -	   	.	(7.157)

Staircase approximation

As already discussed, if the surface interface is z-invariant (entirely or piecewisely), the integration of the system of ordinary differential equations along z can be done via eigenvalue/eigenvector technique, eq.( 7.1 ), because the M-matrix containing the coefficients of the differential equations does not depend on z. The enormous interest in this approach can be explained by the simple technique of integration, much easier to understand

Appendix 7.B: Inverted S-matrix propagation algorithm

In Appendix A we have seen how to avoid numerical instabilities due to the artificially growing terms that appear when the propagation of the field amplitudes from one interface to another is made in the wrong direction, a typical property of a half of the field amplitudes used in the transmission matrix approach.

In some cases (for example, the Integral method applied to multilayer grating, but also coordinate transformation method used for a stack containing different profiles), the numerical solution that has been obtained provides a link, having a form inverse to eq.(7.1 )

Of course, it is easy to obtain the form of eq. ( 7.1 ) by simply inverting p 1 T +  , but better to avoid this, because some blocks of the matrix contain large terms compared to the others. In particular, the block p 1 T ++ +  is responsible for a physical 'antipropagation' from p+1 to p, so that potentially it contains growing terms (as it was that case with p 1 T -- + ) in Appendix A.

We can avoid the direct inversion of p 1 T +  by applying a similar procedure as in Appendix 7.A in order to obtain the S-matrix of the stack. Equation (7.1 ) is expanded in the form:

On the other hand, from eq.( 7.