Fadi Baïda

Abderrahmane Belkhir

Kirill Cherednichenko

Shane Cooper

Richard Craster

Guillaume Demesy

John Desanto

Gérard Granet

Editor Evgeny Popov

Boris Gralak

Leonid I Goray
email: goray@spbau.ru

Sébastien Guenneau

Lifeng Li

Daniel Maystre

André Nicolet

Gunther Schmidt
email: schmidt@wias-berlin.de

Brian Stout

Fréderic Zolla

Theory and Numeric Applications Second Revisited Edition Tryfon Antonakakis

Table of Contents:

12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.9.2 Efficiencies and polarization angles of dielectric sine grating . . . . . . 12.9.3 Efficiencies and polarization angles of metallic echelette grating . . . . 12.9.4 Anomalously absorbing Ag shallow-sine grating in the visible . . . . . 12.9.5 Photonic crystals with Au nanorods in the visible-near-IR . . . . . . . 12.9.6 Lossless photonic crystal with circular rods in the near-and mid-IR . . 12.9.7 Al echelle grating coated by MgF 2 in the VUV . . . . . . . . . . . . . 12.9.8 Au off-plane-grazing-incidence blaze grating in soft x-rays . . . . . . . 12.9.9 W/B 4 C multilayer off-plane-grazing-incidence blaze grating in soft-xrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12.9.10 Flight Mo/Si multilayer rough lamellar grating in the EUV . . . . . . . 12.10Appendix A: Derivation of the recursive algorithm for Separating solver . . . . 12.11Appendix B: Derivation of the recursive algorithm for Penetrating solver . . . . 12.12Appendix C: Derivation of the absorption energy for multilayer gratings . . . . 12.13Appendix D: Derivation of the general connection rule between 2D and 1D gratings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Introduction

This work is part of research that has been pursued by the authors over a long period of time for the purpose of developing accurate and fast numerical algorithms, including the commercial packages PCGrate and DiPoG [12.1, 12.2] designed to model multilayered gratings having mostly one-dimensional periodicity (1D), including roughness, and working in all, including the shortest, optical wavelength ranges at arbitrary optical mounts. The boundary integral equation theory or, briefly, integral method (IM) is presently universally recognized as one of the most developed and flexible approaches to an accurate numerical solution of diffraction grating problems (see, e.g., Ref. 12.3 and Ch. 4 and references therein). Viewed in the historical context, this method was the first to offer a solution to vector problems of light diffraction by optical gratings and to demonstrate remarkable agreement with experimental data. This should be attributed to the high accuracy and good convergence of the method, especially for the TM polarization plane. It does not involve limitations similar to those characteristic of the Coupled-Wave Analysis (CWA), and it provides a better convergence. The disadvantages of this method include its being mathematically complicated, as well as numerous "peculiarities" involved in numerical realization. In particular, quasi-periodic Greens functions and their derivatives appearing as kernels in the integral operators require sophisticated lattice sum techniques to evaluate. Moreover, application of the IM to cases of heterogeneous or anisotropic media meets with difficulties; however, with the volume integral method it is possible to overcome these difficulties easily. Nevertheless, it is on the basis of this theory that all the well-known problems of diffraction by periodic and non-periodic structures in optics and other fields have been solved. In many cases it offers the only possible way to follow up in research. The flexibility and universality inherent in the IM, in particular, enable one rather easily to reduce the problem of radiation of Gaussian waves or of a localized source to that of plane-wave incidence, for which scientists all over the world have a set of numerical solutions. Generalizations of the IM have recently been proposed for arbitrarily profiled 1D multilayer gratings [12.4], randomly-rough x-ray-extreme-ultraviolet (EUV) gratings and mirrors [12.5, 12.6], conical diffraction gratings including materials with a negative permittivity and permeability (metamaterials) [12.7, 12.8], bi-periodic anisotropic structures using a variation formulation [12.9], Fresnel zone plates and diffraction optical elements [12.10, 12.11], and two-dimensional (2D) [12.12, 12.13] and three-dimensional (3D) [12.14] photonic crystals (inclusions) of some geometries, among others.

The IM is so pivotal that one can indicate the few areas where it can be modified and improved to solve particular diffraction problems. By convention they are: (1) physical modelchoice of boundary types, boundary conditions, layer and substrate refractive indices, and radiation conditions; (2) mathematical structure-integral representations using potentials or integral formulas and a multilayer scheme; (3) method of approximation and discretizationdiscretization schemes, choice of basis (trial) and test (weighting) functions, and treatment of coincident points and corners in boundary profile curves; (4) low-level details-calculations and optimization of kernel functions, mesh of discretization (collocation) points, quadrature rules, and solution of linear algebraic systems; (5) implementation enhancements-memory caching, other implementation details. A self-consistent explanation of the existing IMs is beyond the primary scope of the present study. The main purpose of this Chapter is to present a complete description in general operator form of the two IMs applied to 1D multilayer gratings working in conical diffraction mounts and in short waves. Our study also includes the calculus of grating absorption in the explicit form and scattering intensity of randomly-rough gratings using Monte Carlo simulations. For other formal IM treatments and their comparisons, one should rather look to the references of this Chapter as well as to Ch. 4 and to references therein.

Various kinds of electromagnetic features of different nature can exist and be explored in complex grating structures: Bragg and Brewster resonances, Rayleigh anomalies and groove shape features, waveguiding and Fano-type modes, etc. In conical diffraction, the influence of possible types of waves can be mixed. For the purposes of this Chapter, we chose three important types, among many others, of diffraction grating problems to include them in Section 12.9 "Examples of numerical results". They are: bare dielectric or metallic gratings of standard groove shapes working in conical diffraction in the resonance domain; shallow high-conductive or dielectric gratings of various boundary shapes, including closed ones, working in different mounts and supporting polariton-plasmon excitation or Bragg diffraction in the visible-infrared range; bare and multilayer gratings working in grazing-conical or near-normal in-plane diffraction in the soft x-ray-EUV range.

Integral method for one-profile gratings in conical diffraction

The present IM designed for the calculation of the efficiency of bulk and multilayer gratings with arbitrary boundary shapes including micro-and nanoroughness and over an extremely wide wavelength range is considered here in a general operator formalism. In this Section, we consider single-boundary integral equations, which involve boundary integrals of the single and double layer potentials and also the normal and tangential derivatives of single layer potentials. Analytical aspects of boundary integral operators are well represented in publications and, the most relevant of them for the present study, are also in following Sections. The fields are assumed time harmonic. Under these conditions, in classical (in-plane) diffraction the Maxwell system of equations reduces to a single Helmholtz equation; therefore fields are represented in the sequel by scalar functions. They would be vector functions with two components in the case of conical (off-plane) diffraction. In the present Section, we are concerned mostly with conical diffraction, including some notes about metamaterials. Classical diffraction is considered as a particular case with some important details for the implementation.

There exist different ways to transform the diffraction problems under consideration to one-dimensional integral equations over the boundary profile curve of the grating. It is beyond the scope of this chapter to describe the history of applying integral methods to grating problems and the variety of corresponding integral formulations. It should be mentioned that those methods were mostly developed by specialists in physics and optics, and, it seems, they were not aware of the rapid progress in the fields of "boundary integral equations" and "boundary element methods" made in the mathematical community since 1980.

Maxwell equations

We denote by e x , e y and e z the unit vectors of the axis of the Cartesian coordinates. The grating is a cylindrical surface whose generatrices are parallel to the z-axis (see Fig. 12.1) and whose cross section is described by the curve Σ (Fig. 12.2). We suppose that Σ is not self-intersecting and dperiodic in x-direction. The grating surface is the boundary between two regions G ± × R ⊂ R 3 which are filled with materials of constant electric permittivity ε ± and magnetic permeability µ ± . We deal only with time-harmonic fields; consequently, the electric and magnetic fields are represented by the complex vectors E and H, with a time dependence exp(-iωt) taken into account. The wave vector k + of the incident wave in G + × R is in general not perpendicular to the grooves (k + • e z ̸ = 0). Setting k + = (α, -β , γ) the surface is illuminated by an electromagnetic plane wave E i = p e i(αx-β y+γz) , H i = s e i(αx-β y+γz) , (12.1) which due to the periodicity of Σ is scattered into a finite number of plane waves in G + × R and possibly in G -× R. The wave vectors of these outgoing modes lie on the surface of a cone whose axis is parallel to the z-axis. Therefore, one speaks of conical diffraction. The components of k + satisfy β ∈ R and α 2 + β 2 + γ 2 = ω 2 ε + µ + . Note that this condition is satisfied by dielectric media with ε + > 0, µ + > 0 as well as negative index materials, satisfying ε + < 0, µ + < 0. The wave vector k + is expressed using the incidence

angles |θ |, |ϕ | < π/2 (α, -β , γ) = ω √ ε + √ µ + (sin θ cos ϕ , -cos θ cos ϕ , sin ϕ ) .
Note that β > 0 if ε + > 0, µ + > 0, whereas β < 0 for negative index materials. In-plane diffraction corresponds to k + • e z = 0, the case ϕ ̸ = 0 characterizes conical diffraction.

Helmholtz equations

Since the geometry is invariant with respect to any translation parallel to the z-axis, we make the ansatz for the total field (E, H)(x, y, z) = (E, H)(x, y) e iγz (12.2) with the vector functions E, H : R 2 → C 3 . This transforms the time-harmonic Maxwell equations in R 3

∇ ∇ ∇ × E = iω µH and ∇ ∇ ∇ × H = -iωεE , ( 12.3)
with piecewise constant functions ε(x, y) = ε ± , µ(x, y) = µ ± for (x, y) ∈ G ± , into a two-dimensional problem. Indeed, in regions with constant ε and µ we have the equations

∇ γ × E = iω µH , ∇ γ × H = -iωεE , ∇ γ • E = ∇ γ • H = 0 , ( 12.4)
with ∇ γ = (∂ x , ∂ y , iγ). Then from (12.3)

∇ γ × (∇ γ × E) = ω 2 ε µE , ∇ γ × (∇ γ × H) = ω 2 ε µH . (12.5)
Introducing the transverse components

E T = E -E z e z , H T = H -H z e z ,
we derive using (12.4)

∇ γ × (∇ γ × E T ) = γ 2 (E T + E z e z ) + iω µ∇ × (H z e z ) ∇ γ × (∇ γ × H T ) = γ 2 (H T + H z e z ) -iωε∇ × (E z e z )
and

∇ γ × (∇ γ × E z e z ) = (iγ∂ x E z , iγ∂ y E z , -(∂ 2 x + ∂ 2 y )E z ) ∇ γ × (∇ γ × H z e z ) = (iγ∂ x H z , iγ∂ y H z , -(∂ 2 x + ∂ 2 y )H z )
Thus, comparing the components in (12.4) we derive

(ω 2 ε µ -γ 2 )E T = iγ∇E z + iω µ∇ × (H z e z ) ,
(ω 2 ε µγ 2 )H T = iγ∇H z -iωε∇ × (E z e z ) . (12.6) We denote

κ 2 = ε µ - γ 2 ω 2 = ε µ -ε + µ + sin 2 ϕ , ( 12.7)
and conclude from (12.6) that under the condition κ 2 ̸ = 0, which will be assumed throughout, the components E z , H z determine the electromagnetic field (E, H). Furthermore, comparing the third components we derive the relations

ω 2 ε µE z = γ 2 E z -(∂ 2 x + ∂ 2 y )E z , ω 2 ε µH z = γ 2 H z -(∂ 2 x + ∂ 2 y )H z ,
thus E z and H z are solutions of the two-dimensional Helmholtz equations in G ± ∆u + ω 2 κ 2 u = 0 , (12.8) where ∆ = ∂ 2 x + ∂ 2 y denotes the Laplace operator in R 2 . Denote by n = (n x , n y , 0) and t = e z × n = (-n y , n x , 0), respectively, the unit vectors of the normal and the tangent on the surface Γ = Σ × R. Then one finds from (12.6) that

n × E = n × E T -E z t = i ω 2 κ 2 ( γ n × ∇E z + ω µn × (∇ × (H z e z ))
) -E z t

and

n × H = n × H T -H z t = i ω 2 κ 2 ( γ n × ∇H z -ωεn × ∇ × (E z e z ) ) -H z t
Thus the continuity of the tangential components n × E and n × H on the surface Γ = Σ × R leads to the jump conditions for E z , H z across Σ of the form

[ E z ] Σ = [ H z ] Σ = 0 , [ γ ω 2 κ 2 ∂ t H z + ωε ω 2 κ 2 ∂ n E z ] Σ = [ γ ω 2 κ 2 ∂ t E z - ω µ ω 2 κ 2 ∂ n H z ] Σ = 0 . (12.9)
Here ∂ n = n x ∂ x + n y ∂ y and ∂ t = -n y ∂ x + n x ∂ y are the normal and tangential derivatives on Σ, respectively, and [ u] Σ denotes the jump of the function u across the curve Γ. The z-components of the incoming field E i z (x, y) = p z e i(αx-β y) , H i z (x, y) = s z e i(αx-β y) (12.10) are α-quasi-periodic in x of period d, i.e., they satisfy a Floquet condition u(x + d, y) = e idα u(x, y) .

In view of the periodicity of ε and µ this motivates to seek α-quasi-periodic solutions E z , H z . Furthermore, the diffracted fields must remain bounded at infinity, which lead to the outgoing wave condition (OWC). If the profile curve Σ is contained in the strip {(x, y) : |y| < H}, then the quasiperiodicity of solutions implies outside the strip Rayleigh series expansion of the form

(E z , H z )(x, y) = (E i z , H i z ) + ∑ n∈Z (E + n , H + n ) e i(α n x+β + n y) , y ≥ H, (E z , H z )(x, y) = ∑ n∈Z (E - n , H - n ) e i(α n x-β - n y) , y ≤ -H, ( 12.11)
with unknown Rayleigh coefficients E ± n , H ± n ∈ C, and α n , β ± n given by the relations

α n = α + 2πn d , (β ± n ) 2 = ω 2 κ 2 ± -α 2 n .
Then the functions are bounded for |y| → ∞, if we choose the branch of the square root such that Im β ± n ≥ 0, i. e. we set z 1/2 = r 1/2 exp(iφ/2) for z = r exp(iφ), 0 ≤ φ < 2π. In the following it is always assumed that the material in G + satisfies either ε + , µ + > 0 or ε + , µ + < 0, and that the material parameters of the substrate are nonzero complex values with nonnegative imaginary part Im ε -, Im µ -≥ 0. Thus, besides the usual optical materials also interesting negative index materials with ε, µ < 0 are allowed.

If 0 ≤ arg κ 2 ± < 2π, i.e. arg(ε ± + µ ± ) < 2π, then β ± n with 0 ≤ arg β ± n < π are properly defined for all n. However, if arg(ε ± + µ ± ) = 2π, i.e. ε ± , µ ± < 0, we set

κ ± = - ( ε ± µ ± - γ 2 ω 2 ) 1/2 , β ± n = - ( ω 2 κ 2 ± -α 2 n ) 1/2 (12.12)
Summarizing, in case of off-plane diffraction the Maxwell system of equations reduces to two-dimensional Helmholtz equations for vector functions of two components (E z , H z ), which satisfy the OWC (12.11) and are coupled by the jump conditions (12.9).

For in-plane diffraction (γ = 0) we derive from (12.9) the well-know fact, that one can consider the two fundamental cases of polarization separately, i.e. the TE mode (with the z component E z of the electric field E parallel to the grating grooves) and the TM mode (with the z component H z of the magnetic field H parallel to the grating grooves) (see Ch. 2).

Then the surface is illuminated by an electromagnetic plane wave E i = p e i(αx-β y) , H i = s e i(αx-β y) , where α = k + sin θ , β = k + cos θ , with the incidence angle |θ | < π/2, and k + = ω 2 ε + µ + denotes the wavenumber inside G + × R.

The z-components of the total fields E z (TE polarization) or H z (TM polarization) satisfy the Helmholtz equation in G ± except the boundary Σ (12.13) and satisfy the continuity conditions u + Σ = u -Σ = 0, ∂ n u + Σ = q∂ n u -Σ , (12.14) where q = µ -/µ + or q = ε -/ε + for the E z -or H z -component, respectively. In addition, they are subject to OWC (12.11) with β ± n given by relations

(∆ + k 2 ± ) u ± = 0 , k 2 ± = ω 2 ε ± µ ±
β ± n = √ k 2 ± -α 2 n , 0 ≤ arg β ± n < π. (12.15)
Remark 12.2.1 The results of this section can be easily generalized to more general diffraction gratings, where the outer domains G + and G -are bounded by surfaces Σ + and Σ -, surrounding an inner periodic grating structure. For that structure the d-periodic in x material parameters ε(x, y) and µ(x, y) are piecewise continuous functions. Then the z-components E z and H z satisfy Helmholtz equations (12.8) with variable κ 2 = ω 2 ε µ, where ε(x, y) or µ(x, y) are continuous, and meet the jump conditions (12.9) at any discontinuous curve of ε(x, y) or µ(x, y).

Boundary integral operators

Here we describe the application of some mathematical results on boundary integral equations to the solution of the present Helmholtz problems in G ± . Let the common boundary Σ of G - and G + be given by a piecewise C 2 parametrization

σ (t) = (X(t),Y (t)), X(t + 1) = X(t) + d, Y (t + 1) = Y (t) , t ∈ R , ( 12.16)
i.e. the continuous functions X,Y are piecewise C 2 and σ (t

1 ) ̸ = σ (t 2 ) if t 1 ̸ = t 2 .
If the profile Σ has corners, then we suppose additionally that the angles between adjacent tangents at the corners are strictly between 0 and 2π. The potentials which provide α-quasi-periodic solutions of the Helmholtz equation

∆u + k 2 u = 0 with 0 ≤ arg k 2 < 2π (12.17)
are based on the quasi-periodic fundamental solution of period d

Ψ k,α (P) = lim N→∞ i 4 N ∑ n=-N H (1) 0 ( k √ (X -nd) 2 +Y 2
) e iαnd , P = (X,Y ) , (12.18) with the Hankel function of the first kind H

(1) 0 . The series (12.18) converges uniformly over compact sets in

R 2 \ ∪ n∈Z {(nd, 0)} if the condition k 2 ̸ = α 2 n = ( α + 2πn d ) 2 for all n ∈ Z (12.19)
is satisfied. At any point Q = (nd, 0) the fundamental solution has a logarithmic singularity .20) where in the case k < 0 according to (12.7) the fundamental solution is given by (12.20) with

Ψ k,α (P -Q) ≍ e iα(X-nd)-sin(X-nd)) 2π log 1 |P -Q| for P = (X,Y ) near Q. Moreover, setting β n = √ k 2 -α 2 n (recall that Im β n ≥ 0) Poisson's summation formula leads to the representation Ψ k,α (P) = lim N→∞ i 2d N ∑ n=-N e iα n X+iβ n |Y | β n , ( 12
β n = - ( k 2 -α 2 n ) 1/2 .
The single and double layer potentials are defined by

S φ(P) = 2 ∫ Γ φ(Q) Ψ k,α (P -Q) dσ Q , Dφ(P) = 2 ∫ Γ φ(Q) ∂ n(Q) Ψ k,α (P -Q) dσ Q , (12.21)
where Γ is one period of the interface Σ, i.e. Γ = {σ (t) : t ∈ [t 0 ,t 0 + 1]} for some t 0 . In (12.21) dσ Q denotes the integration with respect to the arc length and n(Q) is the normal to Σ at Q ∈ Σ pointing into G -. Obviously, for α-quasiperiodic densities φ on Σ the potentials S φ, Dφ are α-quasiperiodic in X and do not depend on the choice of Γ. They are solutions of the Helmholtz equation (12.17) in G ± and satisfy the radiation condition

u(x, y) = ∞ ∑ n=-∞ u n e iα n x+iβ n |y| , |y| ≥ H . (12.22)
The potentials provide the usual representation formulas. Any α-quasiperiodic solution u of (12.17) in G + satisfying (12.22) admits the representation 1 2

( S ∂ n u -Du ) = { u in G + , 0 in G -, (12.23)
where the normal n points into G -. Under the same assumptions for a function u in G -the representation 1 2

( Du -S ∂ n u ) = { 0 in G + , u in G -, (12.24)
is valid.

Restriction of the potentials S and D to the profile curve Σ are the so called boundary integral operators. The potentials provide the usual jump relations of classical potential theory. The single layer potential is continuous across Σ (S φ) + (P) = (S φ) -(P) = V φ(P) , (12.25) where the upper sign + resp.denotes the limits of the potentials for points in G ± tending in non-tangential direction to P ∈ Σ, and V is a integral operator with logarithmic singularity

V φ(P) = 2 ∫ Γ Ψ k,α (P -Q) φ(Q) dσ Q , P ∈ Σ .
The double layer potential has a jump if crossing Γ:

( Dφ ) + = (K -I)φ, ( Dφ
) -= (K + I)φ (12.26)
with the boundary double layer potential

Kφ(P) = 2 ∫ Γ φ(Q) ∂ n(Q) Ψ k,α (P -Q) dσ Q + (δ (P) -1)φ(P) .
Here δ (P) ∈ (0, 2), P ∈ Σ, denotes the ratio of the angle in G + at P and π, i.e., δ (P) = 1 outside corner points of Σ. The normal derivative of S φ at Σ exists outside corners and has the limits

( ∂ n S φ ) + = (L + I)φ, ( ∂ n S φ ) -= (L -I)φ, (12.27)
where L is the integral operator on Γ with the kernel ∂ n(P) Ψ k,α (P -Q),

Lφ(P) = 2 ∫ Γ φ(Q) ∂ n(P) Ψ k,α (P -Q) dσ Q , P ∈ Σ . (12.28)
Boundary integral methods for second order partial differential equations usually employ also normal derivatives of the double layer potential. The kernel function of this integral operator has a strong singularity of the form |P -Q| -2 near Q = (nd, 0), and must be interpreted as hypersingular or integro-differential operator. Since the computation of this kernel function is rather complicated and time-consuming integral methods for diffraction gratings avoid this operator. However, in the following we need also the tangential derivative of single layer potentials

∂ t ( V φ ) (P) = 2 ∂ t ∫ Γ Ψ k,α (P -Q) φ(Q) dσ Q , P ∈ Σ .
Interchanging differentiation and integration leads to an integral kernel with the non-integrable main singularity t(P)

• (P -Q) |P -Q| 2 ,
where t(P) denotes the tangential vector to Σ at P. Therefore the tangential derivative of single layer potentials cannot be expressed as a usual integral. But it can be interpreted as the Cauchy principal value or singular integral

Jφ(P) = 2 lim δ →0 ∫ Γ\Γ(P,δ ) φ(Q) ∂ t(P) Ψ k,α (P -Q) dσ Q = ∂ t ( V φ ) (P) , ( 12.29)
where Γ(P, δ ) is the subarc of Γ of length 2δ with the mid point P. Similarly, one can define the singular integral

Hφ(P) = 2 lim δ →0 ∫ Γ\Γ(P,δ ) φ(Q) ∂ t(Q) Ψ k,α (P -Q) dσ Q , ( 12.30)
which by using integration by parts gives for α-quasiperiodic φ

Hφ(P) = -2 ∫ Γ Ψ k,α (P -Q) ∂ t φ(Q) dσ Q = -V ( ∂ t φ ) (P) , P ∈ Σ . (12.31) Note that V ∂ t V = V J = -HV .
The integral formulation for the diffraction problems can be derived from the so-called direct or indirect approaches. The direct approach uses the representation formulas (12.23) or (12.24) together with the boundary values (12.25) and (12.26) to derive the boundary integral relations

V ∂ n u + -(I + K)u + = 0 or V ∂ n u -+ (I -K)u -= 0 (12.32)
for quasiperiodic solutions u ± of the Helmholtz equations (12.17) in G ± satisfying (12.22).

Here the unknowns u ± and ∂ n u ± on the profile curve Γ have a direct physical meaning. For the indirect approach, the solution is sought as single or double layer potential with some unknown density.

An important ingredient for discussing the equivalence of integral formulations with the electromagnetic problem is the so-called "Uniqueness Theorem", which looks for conditions on Σ and the wave number k such that the solution of Helmholtz equation (12.17) in G + or G -with zero boundary value on Σ is identical to zero in that domain. The uniqueness of this Dirichlet problem is equivalent to the invertibility of the single layer potential operator V and is guaranteed by two sufficient conditions -Im k 2 > 0; -the profile curve Σ satisfies n y (Q) ≤ 0 for all Q ∈ Σ. Hence, only gratings with overhanging profiles are not covered by the last condition, but it is a quite rare case that the Dirichlet problem with zero boundary data has a nontrivial solution. The only known examples, constructed in Ref. 12.15, are boundaries Σ of very exotic form, which will never appear in practice.

Therefore in the following we will always assume that the "Uniqueness Theorem" is valid. Then any quasiperiodic solution of the Helmholtz equation (12.17) in G + or G -satisfying the OWC (12.22) can be uniquely determined via the representation formulas (12.23) or (12.24), respectively, or can be written as single layer potential V φ with a quasiperiodic density φ, which belongs to some Sobolev-type space of functions given on Γ.

Integral equations for the in-plane case

Let us discuss examples of integral formulations for the in-plane diffraction case. Denote the z-components of the incident wave

u i = { E i z for TE-polarization , H i z for TM-polarization .
Then for bare (one-profile) gratings the problem (12.13) (12.14) (12.11) means that one has to find a solution of

∆u ± + k 2 ± u ± = 0 in G ± , ( 12.33)
satisfying continuity conditions on Σ

u -| Σ = (u + + u i )| Σ , ∂ n (u + + u i ) Σ = q∂ n u -Σ , ( 12.34)
and the outgoing wave condition

u + (x, y) = ∞ ∑ n=-∞ c + n e i(α n x+β + n y) for y ≥ H, u -(x, y) = ∞ ∑ n=-∞ c - n e i(α n x-β - n y)
for y ≤ -H. (12.35) An example of the indirect approach is to look for densities φ + and φ -on Γ such that u + = S + φ + and u -= S -φ - (12.36) satisfy (12.33-12.35), where S ± are the single layer potentials with the fundamental solution Ψ k ± ,α . From (12.34), (12.25) and (12.27) one derives two integral equations on Γ

V + φ + -V -φ -= -u i (L + + I)φ + -q(L --I)φ -= -∂ n u i (12.37)
Here L ± are defined by (12.28) with the kernel ∂ n(P) Ψ k ± ,α (P -Q).

The direct approach uses the relations

V + ∂ n u + -(I + K + )u + = 0 , V -∂ n u -+ (I -K -)u -= 0 , ( 12.38)
following from (12.32), where the double layer potential K ± has the kernel function ∂ n(Q) Ψ k ± ,α (P-Q). The integral equation in Ch. 4.2.3 was derived by putting (12.34) in the second equation in (12.38). Then one gets the system of integral equations

V + ∂ n u + -(I + K + )u + = 0 q -1 V -∂ n u + + (I -K -)u + = -q -1 V -∂ n u i -(I -K -)u i (12.39)
for the unknowns u + and ∂ n u + as functions on Γ.

Another equation system with simpler right-hand side can be obtained by the direct method if one assumes that u i is a solution of ∆u + k 2 + u = 0 in G -and satisfies there (12.35). Hence, one gets additionally to (12.38)

V + ∂ n u i + (I -K + )u i = 0 , leading to the relation V + ∂ n (u + + u i ) -(I + K + )(u + + u i ) = -2u i . (12.40)
As before, (12.34) implies two integral equations, but now for the unknowns u -and ∂ n u -

qV + ∂ n u --(I + K + )u -= -2u i V -∂ n u -+ (I -K -)u -= 0 (12.41)
Note that both the direct and indirect approach lead for TE-and TM-polarization to a system of two linear integral equations with two unknowns. This is the standard approach in the boundary integral method for solving so-called transmission problems. Much effort has been spent in the theoretical and numerical analysis of different integral formulations and approximation methods for their effective solution.

It is popular in grating theory to combine the direct and indirect approaches, which results in a single integral equation for each polarization. This idea goes back to D. Maystre, who already in 1972 proposed this new approach (see Ch. 4). Take for example the representations

u + = S + φ + in G + and u -= 1 2 ( D -u --S -∂ n u - ) in G -.
Then by (12.34) and (12.27)

u -= V + φ + + u i , ∂ n u -= q -1 ( (L + + I)φ + + ∂ n u i )
such that the second relation in (12.38) implies the integral equation

( q -1 V -(I + L + ) + (I -K -)V + ) φ + = - ( q -1 V -∂ n u i + (I -K -)u i ) (12.42)
for one unknown density φ + . Another way is to set

u + = 1 2 ( S + ∂ n u + -D + u + ) in G + and u -= S -φ -in G -.
In this case we derive from (12.34) and (12.27)

u + + u i = V -φ -, ∂ n (u + + u i ) = q(L --I)φ -,
such that (12.40) leads to the single integral equation

( qV + (I -L -) + (I + K + )V -) φ -= 2u i . (12.43)
We see that the combined direct-indirect approach can lead to single integral equations for TE-and TM-problems on one-profile gratings. However, contrary to the pure direct or indirect integral method the equations contain products or compositions of boundary integral operators, which can lead to additional numerical difficulties.

Formulas for Rayleigh coefficients

After having solved one of the integral equation systems (12.37) or (12.41) or one of the single equations (12.42) or (12.43) it is easy to determine the complex amplitudes c ± n of the zcomponents of the plane waves (12.35) reflected and transmitted by the grating-the so-called Rayleigh coefficients. We note that the functions u ± (X,Y ) e -iαX for fixed ±Y ≥ H, respectively, are smooth and d-periodic.

Then c ± n e ±iβ ± n Y is simply the n-th Fourier coefficient of u ± (X,Y ) e -iαX , i.e. c ± n = e ∓iβ ± n Y d ∫ d 0 u ± (X,Y ) e -iαX e -2πnX/d dX = e ∓iβ ± n Y d ∫ d 0 u ± (X,Y ) e -iα n X dX , ±Y ≥ H .
The indirect approach leads to simple formulas. Suppose u ± = S ± φ ± , with known den- sity φ ± . Then for (X,Y ) = P

c ± n = e ∓iβ ± n Y d ∫ d 0 S ± φ ± (P) e -iα n X dX = 2 e ∓iβ ± n Y d ∫ Γ φ ± (Q) dσ Q ∫ d 0 e -iα n X Ψ k ± ,α (P -Q) dX .
It follows from (12.20) immediately, that with Q = (x, y)

∫ d 0 e -iα n X Ψ k ± ,α (P -Q) dX = i 2 e -iα n x+iβ ± n |Y -y| β ± n (12.44) such that c ± n = i dβ ± n ∫ Γ e -iα n x∓iβ ± n y φ ± (Q) dσ Q , where Q = (x, y) . (12.45)
Using the direct approach we find

u ± = ± 1 2 ( S ± φ ± -D ± ψ ± )
with known functions φ ± , ψ ± . Then

c ± n = ± e ∓iβ ± n Y 2d ∫ d 0 ( S ± φ ± (P) -D ± ψ ± (P) ) e -iα n X dX = e ∓iβ ± n Y d ∫ Γ ( φ ± (Q) -ψ ± (Q) ∂ n(Q) ) ∫ d 0 e -iα n X Ψ k ± ,α (P -Q) dX dσ Q .
L. Goray and G. Schmidt: Boundary Integral Equation Methods for Conical Diffraction ...

12.13

Hence from (12.44) we get

c ± n = i 2dβ ± n ∫ Γ ( φ ± (Q) -ψ ± (Q) ∂ n ) e -iα n x∓iβ ± n y dσ Q . (12.46)
Let us apply formulas (12.45), (12.46) to the single integral equation (12.43). The Rayleigh coefficients c - n can be determined from

c - n = i dβ - n ∫ Γ e -iα n x+iβ - n y φ -(Q) dσ Q , ( 12.47)
where Q = (x, y) and φ -is the solution of (12.43). For the reflected waves we get

c + n = i 2dβ + n ∫ Γ ( ( q(L --I)φ --∂ n u i ) - ( V -φ --u i ) ∂ n ) e -iα n x-iβ + n y dσ Q . (12.48)

Integral equations for the off-plane case

Now we describe the integral formulation of the conical diffraction (12.8), (12.49), (12.11) for one-profile gratings. To give the same physical dimension to the functions, we use the vacuum impedance

Z v = (µ v /ε v ) 1/2
, where ε v , µ v denote the vacuum permittivity and permeability, re- spectively, and introduce B z = Z v H z . Noting that γ = ω(ε + µ + ) 1/2 sin ϕ the jump conditions (12.9) are rewritten in the form

[ E z ] Σ = [ B z ] Σ = 0 , [ ε ∂ n E z ε v κ 2 ] Σ = - √ ε + µ + ε v µ v sin ϕ [ ∂ t B z κ 2 ] Σ , [ µ ∂ n B z µ v κ 2 ] Σ = √ ε + µ + ε v µ v sin ϕ [ ∂ t E z κ 2 ] Σ . ( 12.49)
Denoting the z-components of the total fields

E z = { u + + E i z u - , B z = { v + + B i z in G + , v - in G -,
the problem (12.8), (12.49), (12.11) can be written so as to find solutions of

∆u ± + ω 2 κ 2 ± u ± = ∆v ± + ω 2 κ 2 ± v ± = 0 in G ± , κ 2 ± = ε ± µ ± -ε + µ + sin 2 ϕ (12.50)
having on Σ the jumps

u -= u + + E i z , ε -∂ n u - ε v κ 2 - - ε + ∂ n (u + + E i z ) ε v κ 2 + = √ ε + µ + ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∂ t v -, v -= v + + B i z , µ -∂ n v - µ v κ 2 - - µ + ∂ n (v + + B i z ) µ v κ 2 + = - √ ε + µ + ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∂ t u -, (12.51)
and satisfying the OWC

(u + , v + )(x, y) = ∞ ∑ n=-∞ (E + n , B + n ) e i(α n x+β + n y) for y ≥ H, (u -, v -)(x, y) = ∞ ∑ n=-∞ (E - n , B - n ) e i(α n x-β - n y)
for y ≤ -H. (12.52) In order to represent u ± and v ± as layer potentials we assume in what follows that the parameters are such that β ± n = (ω 2 κ 2 ± -α 2 n ) 1/2 ̸ = 0 for all n. Similar to the combined approach of the previous Subsection 12.2.4 resulting in the very simple integral equation (12.43) the solutions in G -are sought as single layer potentials

u -= S -w , v -= S -τ (12.53)
with certain auxiliary densities w, τ, whereas the solutions u + and v + are expressed using (12.23), (12.24)

u + = 1 2 ( S + ∂ n u + -D + u + ) , v + = 1 2 ( S + ∂ n v + -D + v + ) in G + . (12.54)
Here we denote by S ± the single layer potential defined on Γ with the fundamental solution Ψ ωκ ± ,α . Correspondingly D ± is the double layer potential over Γ with the normal derivative of Ψ ωκ ± ,α as a kernel function. As in (12.40) we have

V + ∂ n (u + + E i z ) - ( I + K + ) (u + + E i z ) = 2E i z | Σ , V + ∂ n (v + + B i z ) - ( I + K + ) (v + + B i z ) = 2B i z | Σ , ( 12.55)
where V ± denote the boundary single layer potentials

V ± φ(P) = 2 ∫ Γ φ(Q) Ψ ωκ ± ,α (P -Q) dσ Q , P ∈ Σ ,
and the operators K ± and L ± are defined analogously. Since by Eq. (12.27)

u -| Σ = V -w, ∂ n u -| Σ = (L --I)w, v -| Σ = V -τ, ∂ n v -| Σ = (L --I)τ,
we see from Eqs. (12.55) that the jump conditions (12.51) are valid when the unknowns w, τ satisfy the system of integral equations

ε -κ 2 + ε + κ 2 - V + (L --I)w - ( I + K + ) V -w - √ ε v µ + ε + µ v sin ϕ ( 1 - κ 2 + κ 2 - ) V + ∂ t V -τ = 2E i z , µ -κ 2 + µ + κ 2 - V + (L --I)τ - ( I + K + ) V -τ + √ ε + µ v ε v µ + sin ϕ ( 1 - κ 2 + κ 2 - ) V + ∂ t V -w = 2B i z .
(12.56)

Recall that we suppose κ 2 ± ̸ = 0 and ω 2 κ 2 ± -α 2 n ̸ = 0 for all n. For the analytical and numerical treatment of (12.56), it is advantageous to use the relations

V + ∂ t V -= -H + V -= V + J -
(see the definitions (12.29), (12.30)). Then (12.56) becomes a system of singular integral equations

ε -κ 2 + ε + κ 2 - V + (I -L -)w + (I + K + )V -w - √ ε v µ + ε + µ v sin ϕ ( 1 - κ 2 + κ 2 - ) H + V -τ = -2E i z , µ -κ 2 + µ + κ 2 - V + (I -L -)τ + (I + K + )V -τ + √ ε + µ v ε v µ + sin ϕ ( 1 - κ 2 + κ 2 - ) H + V -w = -2B i z .
(12.57)

for which powerful analytical and numerical methods exist.

If the solution of the system (12.56) is found, then the solution of the conical diffraction problem (12.50)-(12.52) can be determined by the relations

u + = - 1 2 ( S + ( ε -κ 2 + ε + κ 2 - (I -L -)w + √ ε v µ + ε + µ v sin ϕ ( 1 - κ 2 + κ 2 - ) J -τ + ∂ n E i z ) + D + V -(w -E i z )
) ,

v + = - 1 2 ( S + ( µ -κ 2 + µ + κ 2 - (I -L -)τ - √ ε + µ v ε v µ + sin ϕ ( 1 - κ 2 + κ 2 - ) J -w + ∂ n B i z ) + D + V -(τ -B i z )
) ,

u -= S -w , v -= S -τ . (12.58)
A detailed mathematical analysis of the system (12.56) is given in Ref. 12.16. In particular, the following properties have been established:

1. The integral equations are equivalent to the Helmholtz system if the operators V + and V - are invertible.

2. If the profile Σ has no corners, then (12.56) is solvable if ε -+ ε + ̸ = 0 and µ -+ µ + ̸ = 0.

3. If the profile Σ has corners, then (12.56) is solvable if ε -/ε + and µ -/µ + / ∈ [-ρ, -1/ρ] for some ρ > 1, depending on the angles at these corners.

4. The solution of (12.56) is unique if Im ε -≥ 0 and Im µ -≥ 0 with Im(ε -+ µ -) > 0 Remark 12.2.2 The one-boundary solver can be used effectively in multilayer grating problems with separating boundaries, i.e., the maximal y value of a given profile is strictly less than the minimal y value of the next profile above (see Se. 12.5.1). In this case, it is possible to determine the diffracted field of the grating by computing scattering amplitude matrices separately for any profile. For each interface between two different materials, the computation of the scattering amplitude matrices corresponds to solving one-boundary conical diffraction problems with plane waves illuminating the interface from above and below.

Efficiency, absorption, and energy balance

In this part, we give formulas for efficiencies and the absorption of bare gratings under oblique incidences and discuss the energy balance.

Efficiencies in conical diffraction

Diffraction efficiencies or far field patterns for the reflected and transmitted fields can easily be found from the corresponding Rayleigh coefficients of the diffracted outgoing waves. Defining the "energy" as the flux of Poynting's vector

P = Re(E × H)/2 (12.59)
through a normed rectangle parallel to the (x, z)-plane, the ratio of the energies of a reflected or transmitted propagating mode and of the incident wave is defined as the efficiency of that diffracted order. Thus, for a propagating plane wave (E, H) = (p, s) e i(k x x+k y y+k z z) , k = (k x , k y , k z ) with |k| 2 = ω 2 ε µ, the energy is proportional to the y-component of Poynting's vector

P y = 1 2 Re(p z s x -p x s z ) .
Since by (12.6)

p x = - 1 |k| 2 -k 2 z ( k x k z p z + ω µk y s z ) , s x = 1 |k| 2 -k 2 z ( ωεk y p z -k x k z s z )
we obtain

P y = ωk y 2(|k| 2 -k 2 z ) (ε|p z | 2 + µ|s z | 2 ) = ωε v k y 2(|k| 2 -k 2 z ) ( ε ε v |p z | 2 + µ µ v |q z | 2 ) (12.60)
for (E, B) = (p, q) e i(k x x+k y y+k z z) with the modification

B = (µ v /ε v ) 1/2 H.
To find relations between the efficiencies in conical diffraction, let E z , B z be a solution of the partial differential formulation of conical diffraction (12.8), (12.49) and (12.11). The expression of the conservation of energy can be derived from a variational equality for E z and B z in a periodic cell Ω H , which has in x-direction the width d, is bounded by the straight lines {y = ±H} and contains Γ. We multiply Eqs.

(∆ + ω 2 κ 2 ) E z = (∆ + ω 2 κ 2 ) B z = 0 in G ± . respectively with ε ε v κ 2 E z and µ µ v κ 2 B z ,
and apply Green's formula in the subdomains Ω H ∩ G ± . Then by using the quasi-periodicity of E z , B z and the jump relations (12.49), one can derive

∫ Ω H ε ε v ( 1 κ 2 |∇E z | 2 -ω 2 |E z | 2 ) + √ ε + µ + ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∫ Γ ∂ t B z E z - ε + ε v κ 2 + ∫ Γ(H) ∂ n E z E z - ε - ε v κ 2 - ∫ Γ(-H) ∂ n E z E z = 0 , (12.61) ∫ Ω H µ µ v ( 1 κ 2 |∇B z | 2 -ω 2 |B z | 2 ) - √ ε + µ + ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∫ Γ ∂ t E z B z - µ + µ v κ 2 + ∫ Γ(H) ∂ n B z B z - µ - µ v κ 2 - ∫ Γ(-H) ∂ n B z B z = 0 , (12.62)
where Γ(±H) denotes the upper and lower straight boundary of Ω H , respectively, and the normal n on Γ(±H) is directed outward. The outgoing wave condition (12.11) implies

∫ Γ(H) ∂ n E z E z = iβ ( |E + 0 | 2 -|p z | 2 + 2i Im ( E + 0 p z e iβ H ) ) + i ∑ n̸ =0 β + n |E + n | 2 e -2H Im β + n , ∫ Γ(-H) ∂ n E z E z = i ∑ n∈Z β - n |E - n | 2 e -2H Im β - n , ( 12
ε + ε v κ 2 + β |p z | 2 - ε + ε v κ 2 + ∑ β + n >0 β + n |E + n | 2 - ε - ε v κ 2 - ∑ β - n >0 β - n |E - n | 2 = - √ ε + µ + ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) Im ∫ Γ ∂ t B z E z , µ + µ v κ 2 + β |q z | 2 - µ + µ v κ 2 + ∑ β + n >0 β + n |B + n | 2 - µ - µ v κ 2 - ∑ β - n >0 β - n |B - n | 2 = √ ε + µ + ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) Im ∫ Γ ∂ t E z B z , which in view of Im ∫ Γ ∂ t B z E z = Im ∫ Γ ∂ t E z B z leads to β κ 2 + ( ε + ε v |p z | 2 + µ + µ v |q z | 2 ) = ∑ β + n >0 β + n κ 2 + ( ε + ε v |E + n | 2 + µ + µ v |B + n | 2 ) + ∑ β - n >0 β - n κ 2 - ( ε - ε v |E - n | 2 + µ - µ v |B - n | 2
) .

(12.64)

Comparing with (12.60), we see that (12.64) relates the energy of the incident wave, which is proportional to the left side of (12.64) with the energies of the reflected and transmitted modes

β + n κ 2 + ( ε + ε v |E + n | 2 + µ + µ v |B + n | 2 ) and β - n κ 2 - ( ε - ε v |E - n | 2 + µ - µ v |B - n | 2
) , (12.65) respectively. Thus, setting the energy of the incident wave

ε + ε v |p z | 2 + µ + µ v |q z | 2 = 1 , ( 12
.66) from (12.64) we derive for lossless gratings that R + T = 1, where R denotes the sum of reflection order efficiencies

R = ∑ β + n >0 β + n β ( ε + ε v |E + n | 2 + µ + µ v |B + n | 2 ) = ∑ β + n >0 η + n (12.67)
and T is the sum of transmission order efficiencies

T = κ 2 + κ 2 - ∑ β - n >0 β - n β ( ε - ε v |E - n | 2 + µ - µ v |B - n | 2 ) = ∑ β + n >0 η - n . (12.68)
Here, the Rayleigh coefficients E ± n and B ± n can be derived using the formulas presented in Section 12.2.5. For example, (12.45) and (12.58) in G -lead to

E - n = i dβ - n ∫ Γ e -iα n x+iβ - n y w(Q) dσ Q , B - n = i dβ - n ∫ Γ e -iα n x+iβ - n y τ(Q) dσ Q , (12.69) with (x, y) = Q ∈ Γ. Further, the direct integral representation in G + implies E + = 1 2 ( S + φ E -D + ψ E ) , B + = 1 2 ( S + φ B -D + ψ B ) ,
with the known functions (cf. (12.58))

φ E = - ( ε -κ 2 + ε + κ 2 - (I -L -)w + √ ε v µ + ε + µ v sin ϕ ( 1 - κ 2 + κ 2 - ) J -τ + ∂ n E i z ) , ψ E = V -(w -E i z ) . φ B = - ( µ -κ 2 + µ + κ 2 - (I -L -)τ - √ ε + µ v ε v µ + sin ϕ ( 1 - κ 2 + κ 2 - ) J -w + ∂ n B i z ) , ψ B = V -(τ -B i z ) .
Thus, from (12.46)

E + n = i 2dβ + n ∫ Γ ( φ E (Q) -ψ E (Q) ∂ n ) e -iα n x∓iβ ± n y dσ Q , B + n = i 2dβ + n ∫ Γ ( φ B (Q) -ψ B (Q) ∂ n ) e -iα n x∓iβ ± n y dσ Q .
(12.70)

Generalization of energy balance for absorbing bare gratings

One of the most important accuracy criteria based on a single computation is the energy balance that can be generalized in the lossy bulk case described in this Subsection. If the grating is perfectly conducting, then the conservation of energy is expressed by the standard energy criterion R = 1, where R is the sum of the reflection order efficiencies.

If the grating is lossless, Im ε -= 0 and Im µ -= 0, then conservation of energy is expressed by a similar energy criterion (see (12.67) and (12.68))

R + T = 1,
where T is the sum of the transmission order efficiencies.

In a general case, if Im ε -̸ = 0 or Im µ -̸ = 0, then T = 0, R < 1, and the remaining part A of the energy is absorbed in the substrate

A + R = 1.
(12.71)

A is called the absorption coefficient or simply the absorption in the given diffraction problem.

Therefore, an important tool to check the quality of the numerical solution for absorbing gratings is the requirement that the sum of the reflected energy and the absorption energy should be equal to the energy of the incident wave. Besides being physically meaningful, the expression (12.71) is very useful as one of numerical accuracy tests for computational codes and especially important in the x-ray and EUV ranges, and also for plasmonics and metamaterials applications, where absorption plays a predominant role. In the lossy case, one needs an independently calculated quantity A to verify (12.71). For such a quantity, we use the absorption integral defined in Ref. 12.7 and derived bellow.

Absorption for bare gratings

To obtain an expression for the absorption energy we apply Green's formula to E z and B z in the domain Ω H ∩ G + , which gives, since the normal n on Γ is exterior for

Ω H ∩ G + ∫ Ω H ∩G + ( |∇E z | 2 -ω 2 κ 2 + |E z | 2 ) = ∫ Γ(H) ∂ n E z E z + ∫ Γ ∂ n E z E z , ∫ Ω H ∩G + ( |∇B z | 2 -ω 2 κ 2 + |B z | 2 ) = ∫ Γ(H) ∂ n E z E z + ∫ Γ ∂ n E z E z .
(12.72)

The outgoing wave condition (12.11) imply (12.63), such that

Im ∫ Γ(H) ∂ n E z E z = -β |p z | 2 + ∑ β + n ≥0 β + n |E + n | 2 , Im ∫ Γ(H) ∂ n B z B z = -β |q z | 2 + ∑ β + n ≥0 β + n |B + n | 2
Since κ + is real, the imaginary parts on the left of Eqs. 12.72 vanish, and therefore

ε + ε v |p z | 2 + µ + µ v |q z | 2 = ∑ β + n ≥0 β + n β ( ε + ε v |E + n | 2 + µ + µ v |B + n | 2 ) + ε + ε v β Im ∫ Γ ∂ n E z E z + µ + µ v β Im ∫ Γ ∂ n B z B z .
Thus, setting as before in (12.66) the energy of the incident wave

ε 0 ε v |p z | 2 + µ 0 µ v |q z | 2 = 1 ,
the sum of reflection order efficiencies R fulfils

R + ε + ε v β Im ∫ Γ ∂ n E z E z + µ + µ v β Im ∫ Γ ∂ n B z B z = 1 ,
i.e., we derive the conservation of energy for absorbing gratings R + A = 1 with the absorption

A = ε + ε v β Im ∫ Γ ∂ n E z E z + µ + µ v β Im ∫ Γ ∂ n B z B z . (12.73) Note that ∂ n E z = ∂ + n E z and ∂ n B z = ∂ + n B z
are the normal derivatives on Γ of the z-components of the total fields in G + , i.e. the sum of the reflected and the incident fields. Using the jump condition (12.51) the formula for A in the case of conical diffraction can be written in the form

A = κ 2 + β Im ( 1 κ 2 - ( ε - ε v ∫ Γ ∂ - n E z E z + µ - µ v ∫ Γ ∂ - n B z B z + 2 √ ε + µ + ε v µ v sin ϕ Re ∫ Γ E z ∂ t B z
) ) .

In terms of the solution w, τ of the integral equations (12.56) the absorption energy is given by the formula

A = κ 2 + β Im ( 1 κ 2 - ∫ Γ ( ε - ε v (L --I)wV -w + µ - µ v (L --I)τ V -τ ) ) + 2κ 2 + sin ϕ β √ ε + µ + ε v µ v Im 1 κ 2 - Re ∫ Γ V -w J -τ .
(12.74)

Efficiencies and absorption for in-plane diffraction

In the special case of in-plane diffraction (ϕ = 0) these formulas provide for lossless gratings

ε + ε v |p z | 2 + µ + µ v |q z | 2 = ∑ β + n >0 β + n β ( ε + ε v |E + n | 2 + µ + µ v |B + n | 2 ) + ∑ β - n >0 β - n β ( µ + µ - ε + ε v |E - n | 2 + ε + ε - µ + µ v |B - n | 2
) , (12.75) and for absorbing gratings we derive the relation

ε + ε v |p z | 2 + µ + µ v |q z | 2 = ∑ β + n >0 β + n β ( ε + ε v |E + n | 2 + µ + µ v |B + n | 2 ) + 1 β Im ( µ + µ - ε + ε v ∫ Γ ∂ n E z E z + ε + ε - µ + µ v ∫ Γ ∂ n B z B z
) .

(12.76)

We get the well-known expressions for efficiencies η ± n and the heat absorption for TE and TM polarizations

η + n (T E) = β + n β |c + n (T E)| 2 |E i z | 2 , η + n (T M) = β + n β |c + n (T M)| 2 |H i z | 2 , η - n (T E) = µ + µ - β - n β |c - n (T E)| 2 |E i z | 2 , η - n (T M) = ε + ε - β + n β |c - n (T M)| 2 |H i z | 2 , A(T E) = 1 β Im µ + µ - ∫ Γ ∂ n E z E z , A(T M) = 1 β Im ε + ε - ∫ Γ ∂ n H z H z .
12.4 Numerical solution of single-boundary problems Progress in algorithms for numerical solutions of 2D and 3D Helmholtz equations in the last two decades has been nearly comparable with that of computer hardware and experimental nanophotonics. Note that some numerical methods, e.g. differential and CWA, but not integral, inherently suffer from ill-conditionness as wavenumbers or order numbers increase. However, algorithms that possess superior convergence properties are less universal with respect to a scatterer geometry and not well-behaved in the high-frequency range. By these reasons, in the present commercial and non-commercial codes based on the IMs, more classical and robust approaches are mostly used, however with a few possible modifications proposed for low λ /d ratio problems and described in Sec. 12.7.

In practice, the convergence and accuracy of efficiency computation using IMs depend significantly on a proper choice of discretization schemes, quadrature rules, and summation methods for the computation of integral kernels. In order to additionally reduce time (up to an order for some problems) for computation matrices of the above operator equations algorithmic enhancements can be applied to one-boundary solvers. It can be done by using, e.g., cache for exponential functions (plane waves) and cache for kernel functions, as described in the following.

Mathematical results for the integral equations

Here we point out some important mathematical aspects of the integral equations (12.42), (12.43) or (12.57), which are the basis for the described IMs. From the theoretical point of view these equations are almost invertible as operators acting between different Sobolev spaces. If the right-hand side of the integral equation has certain smoothness, say r, i.e. it belongs to a Sobolev space H r , then the solution has smoothness r -1 and belongs to H r-1 . The correct choice of the Sobolev spaces and the parameter r is well understood; it depends on the smoothness of the boundary profile.

Informally speaking, the equations correspond to linear operators of the form

(aI + bK)V +C
where V and K are single and double layer potentials with an arbitrary wavenumber including k = 0, a and b are some constants, and C is compact in the pairs of Sobolev spaces H r-1 → H r . The single layer potential V is always invertible as operator H r-1 → H r . If the profile is smooth, then K is compact in H r → H r , hence the equations generate an operator aV +C 1 , which is invertible for almost all parameters in the pairs of Sobolev spaces for all r. Then any results on single layer potentials equations and their approximate solution can be applied.

If the profile has corners, then V is bounded and invertible as operator H r-1 → H r for 0 ≤ r ≤ 1, K is not compact in H r → H r , but bounded in this range; hence one has to study the above operator in the pairs of Sobolev spaces H r-1 → H r for 0 ≤ r ≤ 1. It is invertible except a discrete number of parameter values and one must apply discretization results for both V and for aI + bK. One serious problem arises that the solution is in general not in L 2 , it has singularities of the form O(ρ -δ ), 0 < δ < 1, where ρ is the distance to the closest edge. Then the solution is not finite at the corner points.

This and some other pure mathematical problems will not be discussed in this Chapter, however they have great influence in practice and should be mentioned here. Concerning the numerical solution of the integral equations (12.42), (12.43) or (12.57), one has to consider at least three important theoretical problems:

1. The discretization of V, K, L, H and their products; 2. The convergence of the chosen numerical method under the condition that the kernel functions are exact or computed with some tolerance;

3. The efficient computation of the kernel functions with given accuracy.

For many types of integral operators quite satisfactory solutions to these problems are known, but not all of them have been applied to the diffraction integrals under consideration.

Approximation of integral equations

Here we describe briefly special Nyström and collocation methods used by the authors for solving the equations (12.42), (12.43) or (12.57). In the described realizations, they are rather simple but robust and universal methods with possibly a small number of kernel computations, because this procedure is rather expensive for the diffraction integrals. The discretization of the products of the integral operators appearing in these equations is done with the separate discretization of the integrals

Ax(t) = υx(t) + ∫ 1 0 K(t, s)x(s)ds .
(12.77)

Here we use the parametrization (12.16). Neglecting the factor e iαX(t) the kernel K(t, s) and the unknown function x(s) are periodic. The case υ = ±1 corresponds to the integrals with K and L, whereas υ = 0 for the integral operators with singular kernels V , H and J. The Nyström discretization of A is based on a quadrature rule of the integral

∫ 1 0 x(s)ds ≈ N ∑ j=1 w j x(t j ) ,
where for periodic functions already the rectangular rule w j = q, t j = jq with q = 1/N provides exponential convergence for smooth functions x. The solution {σ k } N k=1 of the linear system

υσ k + q N ∑ j=1 K(t k ,t j )σ j = y(t k ) , for all t k , k = 1, . . . , N ,
is the Nyström approximate solution of the equation Ax(t) = y(t). If a solution of this discrete problem exists and υ ̸ = 0, the case of second order integral equations, then one gets an approximate solution for all t by

υx N (t) = y(t) -q N ∑ j=1 K(t,t j )σ j and therefore x N (t k ) = σ k . The N × N matrix A N = ∥υδ k j + qK(t k ,t j )∥ N j,k=1
is the Nyström discretization of A. For each element of this matrix only one computation of the kernel is necessary, and often some values can be reused. This discretization is accurate for the diffraction integrals on smooth boundaries. However, for integral operators of the first kind, υ = 0, the simple method is not applicable, since the kernels of V , J and H are singular at the diagonal, i.e. the value of K(t k ,t k ) is not defined. There exist various approaches to apply Nytröm's method also to integral equations of the first kind; one of them is reported in Chapter 4. Another efficient realization was developed in Ref. 12.17 to solve boundary integral equations with the single layer potential of the Helmholtz equation, which can be easily adapted to our situation. Then again, only ∼ N 2 computations of the kernel function are necessary and the accuracy of the modified Nyström discretization is determined by the accuracy of the computations of the kernel functions.

The situation is worse for non-smooth profiles. Then the kernels of all integrals are not differentiable at corner points and the solution of the integral equations is singular. Therefore, the usual Nyström methods with the corners among mesh points make no sense, and there exist several proposals to modify or advance Nyström methods. For example, one can form a modified smooth profile curve excluding small neighborhoods of the corner points and construct the Nyström discretization for the integrals on the modified curve. Another practice is to use graded mesh quadratures, that is, quadratures of various types which become increasingly dense near corner points. Then the discretization of the integrals on profiles with corners via the Nyström method can be made accurate; however, compared with the case of smooth profiles, the resulting matrices are excessively large and possibly very ill-conditioned, such that the method can produce inaccurate results. There are various interesting attempts in the mathematical literature to address these difficulties (see Refs. 12.18-12.20). However, it was found that non-accounting (in respect to the usual approaches) of edge and a few other peculiarities of singular or low-convergent integral equations gives accurate and fast results for small wavelengthand height-to-period ratio diffraction problems (see Sec. 12.7.)

To solve the integral equation Ax = y with a collocation method, one has to choose a set of Ñ approximating functions {φ j } and N collocation points 0 ≤ t 1 < t 2 < . . . < t N < 1. The approximate solution is sought in the form

x N (t) = N ∑ j=1 a j φ j (t)
with the unknown coefficients a j to be determined from the collocation equations

Ax N (t k ) = y(t k ) , for all t k .
Thus the collocation discretization of A is given by the N × N-matrix

A N = ∥Aφ j (t k )∥ N j,k=1 .
The approximating functions should be periodic, therefore trigonometric polynomials, periodic splines (piecewise polynomials) or wavelets are good candidates. For smooth profiles we use

N = 2 Ñ + 1 trigonometric monomials φ j (t) = e 2πi jt , j = -Ñ, -Ñ + 1, . . . , Ñ -1, Ñ
for which the expressions Aφ j are cheap to compute. For example, the integrals with singular kernel can be written as

V φ(t) = - ∫ 1 0 log(4 sin 2 π(t -s))φ(s) ds + ∫ 1 0 g V (t, s)φ(s) ds Hφ(t) = ∫ 1 0 cot π(t -s)φ(s) ds + ∫ 1 0 g H (t, s)φ(s) ds ,
where the functions g V , g H are differentiable and periodic in t and s. Thus, the integrals with these kernels can be accurately approximated with the Nyström discretization. For the main parts, the relations

- ∫ 1 0 log(4 sin 2 π(t -s))φ j (s) ds = { φ j (t)/| j| , j ̸ = 0 , 0 , j = 0 , ∫ 1 0 cot π(t -s)φ j (s) ds = { i sign( j)φ j (t) , j ̸ = 0 , 0 , j = 0 , (12.78)
are used, resulting in only one computation of the kernel functions per element of the final collocation matrix. This approach can also be used for profiles with corners, but the accuracy of the trigonometric collocation deteriorates. This is caused by the poor Fourier series approximation of functions, which are singular at corner points. In this case, collocation with splines is advantages, which are able to approximate singular functions on graded meshes. The drawback of splines collocation is that similar to the Nyström discretization graded mesh quadratures are needed to determine accurate collocation discretizations of the integrals. In Sec. 12.4.4 we describe a combination of the trigonometric and spline collocation.

Nyström discretization with modifications

We use the piecewise constant Nyström method (see Subsection before) with the matrix elements such as

A N = ∥υδ ki + qK(t k ,t i )∥ N i,k=1
(12.79)

The present Nyström method proceeds by approximating the values of the current density function φ -of (12.43) at the quadrature nodes (X(t i ),Y (t i )), i = 1, N; by solving the system

φ -(X(t i ),Y (t i )) + N ∑ k=1 c ik φ -(X(t k ),Y (t k )) = b(X(t i ),Y (t i )) , (X(t i ),Y (t i )) ∈ Γ . (12.80) of N linear equations in the N unknowns φ -(X(t i ),Y (t i ))
with composed coefficients c ik containing matrix elements of different operators.

For relatively shallow profiles, the nodes can be uniformly put along the x-coordinates. But the approximately uniform distribution with respect to the arc-length using the parametrization (12.16) is more universal and makes it possible to treat, for example, the lamellar or any other boundary profile with abrupt slopes or non-functions by the integral method without any additional effort on the user side. The principal parameter, with respect to which the convergence is evaluated, is the number N of discretization points on each boundary. In some codes N may vary from one boundary to another, which can be useful. In the present study, let us discard this option for simplicity. Quadratures for operators with continued kernels in our codes are performed by the trapezium (= rectangle) integration rule.

The present numerical solution of the integral equations is based on a simple modification of the Nyström discretization with piecewise constant weighting functions. The choice of a discretization of integral equations such (12.80), PCGrate software default option for most cases, requires a standard regularization of integrals. For smooth curves, the convergence order is determined from the accuracy of computing the fundamental solution, which is N -3 if only the first derivatives of the parametrization of the curve are used. The problem becomes harder if the curve contains corners. Then the double layer potentials K and also potentials L are not compact (the single layer potential is always compact in L 2 ). So the usual Nyström method is very problematic to treat corners; however some interesting modifications can be found in Ref. 12.20, though they are only applicable to K and L. The treatment of the integral kernels K(t k ,t i ) in (12.79) is connected with infinity for the single layer potential in the coincide points, well behaved for the double layer potential on smooth curves, and discontinuous at corner points for the double layer potential. We use three types of single-term corrections: for the single layer potential-taking into account their logarithmic singularities; for the double layer potential-accounting for the profile curvature; and for various kernel functions-acceleration terms applied to the truncated series. So, generally we do modifications in both diagonal and non-diagonal matrix elements in respect to the standard Nyström matrices having N × N regular coefficients.

The integral operators with a weak singularity in the diagonal terms can be split into a sum

V ii = d -1 ∫ Γ ln |2 sin[(P -Q)/2]|φ -dσ Q +CO ii , (12.81)
where CO is a compact operator with continuous kernel. The first integral operator in (12.81) can be calculated easily for some approximations of the current density function ϕ -. It is worth noting that a regularization can be used even at corner nodes of a non-smooth boundary. In the presence of a profile with corners (piecewise linear), the sampling and quadrature nodes are set in such a way that all corners are nodes and the curvature corrections are applied by adding the corner term to the diagonal values. For calculations of shallow gratings having a lot of uniformly-distributed edges (multi-polygonal), another version of the quadrature formula can be applied: the nodes are set in such a way that every corner lies half-way between the nodes adjacent to it and no curvature-like single-term corrections are added (see Sec. 12.7). However, for calculations of deep grating having abrupt and/or long slopes such a simplified approach does not work and a formula involving the normal derivative of the Green function should be used.

The diagonal element of the discretization matrix corresponding to a corner point takes the value (cf. Ref. 12.21, p. 120):

e ii = (2N) -1 (K L (t i ) + K R (t i )) + 1/2 -ζ i /2π , (12.82)
where ζ i is the exterior angle between adjacent tangents at the corner point P = (X(t i ),Y (t i )) and (cf. Ch. 4, Eq. (4.92))

K L,R (t i ) = y ′ (t i ) [ (2d) -1 ∞ ∑ n=-∞ α n /β n + iα 0 /2π ] -y ′′ (t i )[(1 + (y ′ (t i )) 2 /4π] , (12.83)
where K L,R (t i )) means the left-and right-sided limits, respectively, of the kernel function value at the corner point. For gratings with smooth boundaries, these corrections yield the overall error estimate O(N -3 ) for diffraction amplitudes in both polarizations. However, the above simple singularity accounting is insufficient to match such accuracy of the discretization near the corners and a truncation rule together with some acceleration technique should be applied to the truncated kernel function series. Therefore, in computations of kernels, we use a direct summation approach with possible single-term corrections of corresponding matrix elements (see Sec. 12.4.5.1). The matrices of the discretized operators contain the values of the corresponding kernel functions divided by the number N of segments between collocation points. We use in our codes a few different algorithms for solving linear systems of algebraic equations. It can be either the direct Gauss-Jordan elimination method (Gauss) or the non-direct Full orthogonalization method (FOM) which is similar to the Generalized minimum residual method (GMRES). For the FOM case, the number of iterations until a prescribed residual error is reached depends, of course, on the refraction indices, the number of accounting diffraction modes, and the profile shape, but it is nearly independent of the number of unknowns. Note that discretization of the multi-boundary integral equations can be treated by the same simple manner without modifications.

Hybrid trigonometric-spline collocation

Here we describe collocation methods for solving the integral equations (12.57) of conical diffraction, which contain the singular integral H + . The collocation discretization requires the computation of this integral for special basis functions, which is simpler than the Nyström discretization, which uses point values of the strongly singular kernel function of H + .

We consider a parametrization of Γ given by (12.16). In the case of a smooth profile Σ a trigonometric collocation method is used, i.e. we use approximations (12.84) where the coefficients {a k }, {b k } are such that the system (12.56) is satisfied at the

w(σ (t)) e -iαX(t) |σ ′ (t)| ≈ w N (t) = Ñ ∑ k=- Ñ a k e 2πikt , τ(σ (t)) e -iαX(t) |σ ′ (t)| ≈ τ N (t) = Ñ ∑ k=- Ñ b k e 2πikt ,
N = 2 Ñ + 1 collocation points t k = k/N, k = 1, . . . , N.
The advantage of using trigonometric methods is that the integral operators V ± and H + with singular kernels can be approximated properly. For example, using the parametrization σ (t) the single layer potential operator of w can be approximated by

V ± w(σ (t)) ≈ -e iαX(t) ( ∫ 1 0 log(4 sin 2 π(t -s)) w N (s) ds + ∫ 1 0 g ± (t, s)w N (s) ds ) ,
and the singular integral J ± w by

H + w(σ (t)) ≈ -e iαX(t) ( ∫ 1 0 cot π(t -s) w N (s) ds + ∫ 1 0 j + (t, s)w N (s) ds ) ,
where the functions g ± (t, s), j + (t, s) are differentiable and periodic in t and s. The action of the integral operators with the kernels log(4 sin 2 π(ts)) and cot π(ts) on trigonometric polynomials is given analytically; compare (12.78). All other integrals have differentiable kernels and they are approximated by the trapezoidal rule like in the Nyström method described above. So the discretization error depends only on the error made in computing the functions g ± (t, s), j + (t, s) and the continuous kernels of K + and L -, i.e. in computing the fundamental solution and their derivatives. Here we use the exact Ewald method (cf. Section 12.4.5.2) with a number of summation terms to ensure discretization errors of order N -3 . Finally, the operator products V + L -, K + V -or H + V -are approximated by the products of the corresponding discretization matrices.

If the profile curve has corners, then the convergence properties of methods with only trigonometric trial functions deteriorate due to singularities of the densities w and τ of the form O(ρ -δ ), 0 < δ < 1, where ρ is the distance to the closest edge. In boundary element methods it is common to use piecewise polynomial trial functions on meshes graded towards corner points. But due to the complicated form of their kernels the quadrature of the integral operators acting on piecewise polynomials is very expensive. Therefore we use a modification of the trigonometric collocation scheme with a fixed number of piecewise polynomial trial functions.

In the beginning, we introduce meshes of collocation points which contain the corners and are graded towards the corner points. This can be derived by changing the parametrization (12.16), for example, if σ (t j ) is a corner point, then σ ′ (t j ) = σ ′′ (t j ) = 0 implies grading towards the corner. Further, for each collocation point t k there exists a Lagrangian trigonometric polynomial p k (t) of degree N such that

p k (t j ) = δ k j , k, j = 1, . . . , N ,
where δ k j is Kronecker's delta. For each edge and a fixed number of collocation points t k around it we replace the corresponding Lagrangian trigonometric polynomial p k (t) by a cubic spline s k (t) on the graded mesh with s k (t j ) = δ k j . Thus we get a hybrid trigonometric-spline collocation method, which combines the efficient computation of the integrals for trigonometric polynomials with the good approximation properties of piecewise polynomials on graded meshes near edges. The values at the collocation point t j of the integrals on the basis spline s k are computed by a composite Gauss-quadrature with a quadrature mesh geometrically graded towards t j and depending on the distance |σ (t k ) -σ (t j )|. This leads to a fixed number of additional calculations of the fundamental solutions Ψ k ± ,α for each discretization level compared with the pure trigonometric method, which is however compensated by a significant higher accuracy.

Evaluations of kernels

In spite of many research efforts (see, e.g., Refs. 12.22-12.24)-computation of the kernels remains a most time-critical part of integral method for periodic structures. Convergence of the kernels deteriorates significantly as the distance between function's arguments (a discretization point or/and a quadrature node) tends to zero, and especially near edges and at high frequencies. For more discussions we refer the reader to Ch. 4. Some "crash test" calculations on PCGrate codes can be found in Ref. 12.25 and also in numerical examples of Sec. 12.9. The Ewald sum method is quite intricate and widely used (see, e.g., Ch. 6). It is based on a separation of the infinite sum into slowly and rapidly convergent parts and, then, a transformation of the slowly convergent part using the Poisson formula and error functions. Although Ewald methods are proven to be quite efficient for many diffraction grating problems, it has turned out to exhibit poor numerical properties in short waves.

Direct kernel summation

In the described discretization method, we use a direct approach for the evaluation of kernel functions based on Poisson's summation formula (see (12.20)) and a special rule for various kernels with positive and negative summation index. In the simplest case, the series is truncated symmetrically at the lower summation index -P and upper index P; where P is an integer defined by P ≈ gN. (12.85)

For many grating efficiency problems, the number of terms P with plus or minus sign you choose should be fifty percent of a number of discretization points N ("the golden rule" and default value of PCGrate codes; for more see Sec. 12.5.7). In difficult cases like those of highly conducting blazed or very deep gratings, echelles, grazing incidence, and, especially, for the TM polarization one can try to optimize convergence and accuracy by varying g at a given number of discretization points. This parameter can be optimized at small values of N and the ratio is kept constant as N increases. Fortunately, it should be done only in very exceptional cases.

As an easy remedy to accelerate convergence of the series representing the kernels, we use the Aitken δ 2 method [12.26], which is a simple one-term improvement over a popular ac- celeration technique described in Ch. 4. The precision of Aitken's method for individual values of the kernels, especially at close arguments, is inferior to that provided by the Kummer acceleration used in the IESMP code ([12.27]) or by the Euler-Knopp method [12.28], but one would not benefit from extra accuracy in the end. Such more accurate acceleration techniques make sense in combination with higher-order collocation or Galerkin methods. However, it is usually a difficult task to achieve an acceptable combination of robustness and wide-range applicability with higher-order codes. In the Aitken method suppose series S = Σa k has approximately geometric convergence, then the sum

SK+1 = ∑ k≤K a k + a 2 K /(a K -a K+1 ) , (12.86)
is a better converged series than S K+1 having the same number of terms. Both versions of discretization near corners together with the acceleration technique described above are found to yield approximately the same convergence rate O(N -2+ε ), where 0 < ε < 0.5 apparently depends on boundary profile geometries. Such a regularization of the weakly singular and singular integral operators, together with an acceleration of the truncated kernel function series, theoretically and numerically leads to higher rates of convergence and to bounded condition numbers of the discretization matrices. Though, for discretization numbers N of practical interest, no advantage of regularization is observed in our numerical experiments at very small λ /d ratios (see Sec. 12.7), even for smooth boundaries.

Ewald's method

It has turned out that acceleration techniques for the summation approach is not efficient if the second argument y has small modulus |y| (cf. Ref. 12.22), which frequently occur in the quadrature of integrals for graded meshes near corners or for thin layers, i.e., very close profile curves. In this case one can use the following summation algorithm for the integral kernel which is based on Ewald's method (cf. Ref. 12.23), providing a good overview on various methods for the computation of the fundamental solution.

The idea is to split the slowly converging series (12.20) into two quickly converging series. To simplify of presentation we consider the infinite series

Ψ(x, y) = i 4π ∑ n∈Z e inx+iβ n |y| β n (12.87) with β n := √ k 2 -α 2 n and α n := n + α and let Re β n , Im β n ≥ 0. Ψ(x, y) is 2π-periodic in x. Ewald's method is based on the relation ie iβ n |y| β n = ∫ a 2 0 exp ( β 2 n t - y 2 4t ) dt √ πt + i 2β n ( e -iyβ n erfc ( -iaβ n + y 2a ) + e iyβ n erfc ( -iaβ n - y 2a
) )

, which is valid for any a > 0 and β n ̸ = 0. Here

erfc(z) := 2 √ π ∫ ∞ z e -t 2 dt
is the complementary error function. Thus we have Ψ = Ψ e + Ψ w with the two sums

Ψ e (x, y) = 1 4π ∑ n∈Z e inx ∫ a 2 0 e β 2 n t-y 2 /4t dt √ πt , ( 12
.88)

Ψ w (x, y) = i 8π ∑ n∈Z e inx β n ( e -iyβ n erfc ( -iaβ n + y 2a ) + e iyβ n erfc ( -iaβ n - y 2a ) ) . (12.89) Since β 2 n = k 2 -α 2
n , the first sum (12.88) takes the form Denoting r 2 m := (x -2πm) 2 + y 2 and using the series expansion of e k 2 t yields

Ψ e (x, y) = 1 4π ∑ n∈Z e inx ∫ a 2 0 e (k 2 -α 2 n )t-y 2 /4t dt √ πt = 1 4π ∫ a 2 0 e k 2 t
∫ a 2 0 e k 2 t e -r 2 m /4t dt t = ∞ ∑ j=0 k 2 j j! ∫ a 2 0 t j-1 e -r 2 m /4t dt = ∞ ∑ j=0 (ak) 2 j j! E j+1 ( r 2 m 4a 2
) with the exponential integral function E j of degree j

E j (z) := ∫ ∞ 1 e -zt t j dt .
Thus we obtain the representation

Ψ e (x, y) = e -iαx 4π ∑ m∈Z e 2πimα ∞ ∑ j=0 (ak) 2 j j! E j+1 ( r 2 m 4a 2 ) . (12.91) Since E j+1 (z) ≤ e -z z + j , z > 0 ,
the expression (12.91) for Ψ e is quickly converging if (x, y) ̸ = (2πm, 0) with a speed of convergence increasing as the parameter a gets smaller. The function Ψ w can be transformed to a computationally suitable form by using the scaled complementary error function

w(z) := e -z 2 erfc(-iz) = e -z 2 2 √ π ∫ ∞ -iz e -t 2 dt = 2 √ π ∫ ∞ 0 e -t 2 e 2izt dt , (12.92) which has the properties w(-z) = w(z) , w(-z) = 2e -z 2 -w(z) , | w(z)| ≤ 1 for Im z ≥ 0 . (12.93) Using e ∓iyβ n erfc ( -iaβ n ± y 2a ) = e a 2 k 2 e -a 2 α 2 n e -y 2 /4a 2 w ( aβ n ± i y 2a
) ,

we can write (12.89) in the form

Ψ w (x, y) = i e -y 2 /4a 2 e a 2 k 2 8π ∑ n∈Z e inx e -a 2 α 2 n β n ( w ( aβ n + i y 2a ) + w ( aβ n -i y 2a
) )

. ( ) )

= i 4π e i|y|β n β n .
Introducing the finite set P := {n ∈ Z : Im β n + Re β n < |y|/[2a 2 ]}, the function Ψ w is decom- posed into an exponentially converging series and two finite sums

Ψ w (x, y) = i e -y 2 /4a 2 e a 2 k 2 8π { ∑ n∈Z\P e inx e -a 2 α 2 n β n ( w ( aβ n + i y 2a ) + w ( aβ n -i y 2a
) )

+ ∑ n∈P e inx e -a 2 α 2 n β n ( w ( aβ n + i |y| 2a ) -w ( -aβ n + i |y| 2a ) )} + i 4π ∑ n∈P e inx e i|y|β n β n .
(12.96)

In particular, in the case y = 0 which occurs frequently for binary gratings, we obtain the exponentially converging series

Ψ w (x, 0) = i e a 2 k 2 4π ∑ n∈Z e inx e -a 2 α 2 n β n w ( aβ n ) .
Note that the speed of convergence of the series in (12.96) is increasing as the parameter a gets larger.

The representation Ψ = Ψ e + Ψ w is also used for the computation of the gradient of Ψ

(∂ x + iα)Ψ(x, y) = - 1 4π ∑ n∈Z α n e inx+iβ n |y| β n , ∂ y Ψ(x, y) = - 1 4π ∑ n∈Z sign(y)e inx+iβ n |y| ,
which is needed to compute the kernels of the operators K, L, J and

H. Since ∂ z E j (z) = -E j-1 (z) with E 0 (z) := e -z /z, the derivatives of Ψ e are (∂ x + iα)Ψ e (x, y) = - e -iαx 2π ∑ m∈Z (x -2πm) e 2πimα ( e -r 2 m /4a 2 r 2 m + ∞ ∑ j=1 (ak) 2 j 4a 2 j! E j ( r 2 m 4a 2
) ) ,

∂ y Ψ e (x, y) = - y e -iαx 2π ∑ m∈Z e 2πimα ( e -r 2 m /4a 2 r 2 m + ∞ ∑ j=1 (ak) 2 j 4a 2 j! E j ( r 2 m 4a 2
) ) .

(12.97)

The derivatives of Ψ w are given by

(∂ x + iα)Ψ w (x, y) = - e -y 2 /4a 2 e a 2 k 2 8π { ∑ n∈Z\P α n e inx e -a 2 α 2 n β n ( w ( aβ n + i y 2a ) + w ( aβ n -i y 2a
) ) ) ) )

+ ∑ n∈P α n e inx e -a 2 α 2 n β n ( w ( aβ n + i |y| 2a ) -w ( -aβ n + i |y| 2a ) ) } - 1 4π ∑ n∈P α n e inx e
= iβ n e -y 2 /4a 2 ( w

( aβ n -i y 2a ) -w ( aβ n + i y 2a
) ) .

The numerical calculation of the exponential integral E j and its derivatives and of the scaled complementary error function w(z) present no problem using standard routines. The value of the parameter a should be chosen small enough to ensure the rapid convergence of the series for Ψ e and its derivatives and large enough to ensure the rapid convergence of the series representations for Ψ w and its derivatives. After numerical tests we found that the choice a|k| = 6 is a good compromise.

Cache for exponential functions (plane waves)

In order to reduce time for computation matrices of the above operator equations, a simple but effective acceleration was implemented in authors' codes at the algorithmic level, i.e. cache for exponential functions (plane waves). It assumes a large time-memory tradeoff. The amount of memory required for cache can be calculated in advance in each case and adjustments (cache off or partial) are done automatically. More acceleration can be reached in some cases, e.g. if one uses cache for kernel functions (see Sec. 12.6.2). Calculation of kernel functions makes extensive use of typical multiplicative combinations of exponential functions

exp{iα n (X i -X k ) + iβ n |Y i -Y k |} =    E - n,i /E + n,k if Y i ≥ Y k ; E + n,i /E - n,k if Y i < Y k .
(12.100)

Here E ± k,i = exp{iα n X i ± iβ n Y i }. (12.101)
Let N be the number of discretization points on a given boundary, that is, the subscripts i and k in the above expressions assume N different values. Let P be the number of exponential terms to be stored in the cache. That is, the index n assumes P values situated symmetrically (with a possible ±1 imbalance) with respect to 0 (see (12.85)). Normally P is the maximum number of negative or positive exponential terms used in computations of kernel functions. If, however, there is not enough fast memory in the system, a partial cache is used, where some exponents are pre-computed and extracted form cache in the course of the kernel function computations, while other exponents are evaluated on the spot. In total, 2 PN exponents are pre-computed for every boundary. The value of P may vary, depending on which (if at all) acceleration method is used for the series summation for the kernel functions, however for many cases P ≤ N . So, memory expenditure is again of the order N 2 per layer. The pre-computed exponents share the same memory for every layer, so newer values override old ones. Unlike the kernel function cache (see Sec. 12.6.2), saving the precomputed exponents for a potential re-use in further layers with same refractive indexes does not make much sense: pre-computation only needs O(N 2 ) operations per layer, which is a tiny fraction of the total, which is of the order N 3 .

Keeping track of the stored elements order this case does not call for any special technique like as binary trees: a two-dimensional array is all one needs. However, a difficulty of another sort pops up. The numbers β n have nonzero imaginary parts when |n| exceeds some n 0 , and the asymptotics of Im β n is linear as |n| grows indefinitely. Depending on the signs of Y i in Eq. (12.101), the exponents easily go beyond the underflow and overflow limits in the standard floating-point arithmetic. However, the absolute values of resulting ratios (12.100) are always not greater than 1.

To resolve this problem, the data {E k,i } are stored in the format {mantissa, order}; see Ref. 12.29. The order is represented by a variable of an integer type. It can be unusually large (positive or negative) if one thinks about typical orders in engineering calculations. For example, the values Im β n = 1000 and Y i = 10, though rather extreme, can occur in grating calculations. But for the data structure we describe, numbers like exp (10 4 ) are nothing unusual and totally within its capacity. In our program, the 2-byte C type short int is chosen for orders, which suffices for all practical purposes.

We fix a huge positive B (the "base"); in the program B = 10 20 , a more or less arbitrary value. Every nonzero real or complex number Q is then uniquely represented in the form

Q = B q • M, 1 ≤ | M| < B (12.102)
with integer q. The only arithmetical operation needed for (12.100

) is division Q/Q ′ given that |Q| ≤ |Q ′ |. Assuming Q ′ = B q ′ • M ′ , set Q Q ′ =      M/M ′ , if q = q ′ , ( M/M ′ )B -1 , if q = q ′ -1, 0, if q < q ′ -1. (12.103)
The divisions on the right are carried out in the in standard floating-point format.

Solving diffraction of multilayer gratings

The use of coatings has many applications in diffraction gratings. We shall consider two algorithms for conical diffraction by multilayer gratings based on the integral methods for oneprofile gratings, which are theoretically able to deal with multilayered gratings without limitations concerning the shape of the interfaces or the conductivity of the layers. The choice of a numerical method to solve the multi-boundary integral equations is to a large extent independent of other implementation details of the single-boundary algorithm. It is not even necessary to use the same method for every boundary, provided that adjacent boundary solvers have a common data interface. In the hope of making the algorithm more accessible, we explicitly write out a chain of operator equations to emphasize the upper-level structure of the multilayer algorithms. Details which are not pertinent to the structural level are omitted here, but are well discussed in other Sections and Appendices. Assuming the potential operators are available as ready-to-use building blocks, an object-oriented implementation of the operator equations becomes relatively easy. Our description of the multilayer schemes below emphasizes its structural aspects from the perspective of an object-oriented implementation. There are two different multilayer solvers implemented in the authors' codes: the 'Separating' multi-boundary solver based on the scattering amplitude matrix algorithm described in Appendix A and the 'Penetrating' multi-boundary solver based on recursive marching algorithms described in Appendix B. The first one is restricted to multilayer gratings with horizontally separated boundary profiles, where it is possible to define a plane layer in between that does not cross the upper or the lower interface. Then one can use plane-wave Rayleigh expansions of the electromagnetic field between the interfaces and work with the scattering matrices for that interface. The second algorithm works in the case of interpenetations of interfaces, but is numerically more expensive than the first one.

Mathematical aspects of multi-boundary integral operators are nontrivial, however well represented in this Chapter and many publications. For example, transparent and detailed exposition, including a discussion of various marching schemes that avoid hypersingular potential operators, is given in Ref. 12.21.

Gratings with separating boundaries

Let us now consider a multilayer diffraction grating with period d formed by a stack of M relief and/or rod gratings characterized by grating profiles Σ j , j = 0, . . . , M -1. More precisely, the structure consists of material layers which are separated by continuous profiles and may contain rod gratings. The different media are numbered from top to bottom; see Figure 12.3, G 0 and G M are the semi-infinite top and bottom layers. To apply a scattering matrix approach, we assume that the interfaces Σ 0 , . . . , Σ M-1 between the M + 1 homogeneous material domains G 0 , . . . , G M are separated, i.e. between adjacent interfaces Σ j and Σ j-1 there exists a strip {u j < y < d j-1 } not crossing the interfaces. The structure of the multi-profile grating is characterized by the permittivity and permeability functions ε(x, y) and µ(x, y), which are constant on the domains G j . Its values in G 0 and G M are denoted by ε 0 , ε M and µ 0 , µ M , respectively. Further we denote

G M G 2 G 2 G 2 G 1 G 0 G 3 y x . . . . . d Σ Σ Σ 3 Σ 0 Σ Σ 2 3 Σ M-1 1 1 1 u d 2
κ 2 0 = ε 0 µ 0 cos 2 ϕ , κ 2 M = ε M µ M -ε 0 µ 0 sin 2 ϕ .
As in the case of one interface, the z-components E z ,

B z = (µ v /ε v ) 1/2 H z satisfy Helmholtz equations ( ∆ + ω 2 κ 2 ) E z = ( ∆ + ω 2 κ 2 ) B z = 0 (12.104)
κ 2 = ε µε 0 µ 0 sin 2 ϕ , in the domains G j and the transmission conditions at the interfaces Σ j

[ E z ] Σ = [ B z ] Σ = 0 , [ ε ∂ n E z ε v κ 2 ] Σ = - √ ε 0 µ 0 ε v µ v sin ϕ [ ∂ t B z κ 2 ] Σ , [ µ ∂ n B z µ v κ 2 ] Σ = √ ε 0 µ 0 ε v µ v sin ϕ [ ∂ t E z κ 2 ] Σ . (12.105)
The light is incident from G 0 and we are interested in the Rayleigh coefficients E ± n , B ± n of the series expansions

(E z , B z )(x, y) = (E i z , B i z ) + ∑ n∈Z (E + n , B + n ) e i(α n x+β (0) n y) , y ≥ H, (E z , B z )(x, y) = ∑ n∈Z (E - n , B - n ) e i(α n x-β (M) n y) , y ≤ -H, ( 12.106)
where the half spaces {y ≥ H} and {y ≤ -H} are contained in the semi-infinite layers G 0 and G M , respectively. According to (12.7) we have β

( j) n = (ω 2 κ 2 j -α 2 n ) 1/2 if 0 ≤ arg(ε j + µ j ) < 2π and β ( j) n = - ( ω 2 κ 2 j -α 2 n ) 1/2 if ε j , µ j < 0.
We study the off-plane diffraction for gratings with separated interfaces using the robust algorithm (for the derivation, see App. A) for modeling layered gratings (an overview is given, for example, in Ref. 12.30). The present method extends the S-matrix algorithm given by D. Maystre in Ref. 12.31 for the integral method and in-plane diffraction. It is a recursive algorithm to determine operators R 0 and T 0 , which map the coefficients (p z , q z ) of the incoming plane wave (E i z , B i z ) = (p z , q z ) e i(αx-β y) to the vectors of Rayleigh coefficients {(E + n , B + n )} n∈Z of the reflected and {(E - n , B - n )} n∈Z of the transmitted fields, cf. (12.11). To this end, the multi-profile problem is split into simpler scattering problems for one-profile gratings, which are formed by the profiles Σ j and separate optical materials with the parameters ε j , µ j and ε j+1 , µ j+1 .

We give a formal operator description of the marching procedure for R 0 and T 0 . For each profile there exist scattering operators, which map the Rayleigh coefficients of an incoming field to the Rayleigh coefficients of the reflected and transmitted fields. More precisely, the grating with profile Σ j diffracts the α-quasi-periodic incoming field

∑ n∈Z (A j n ,C j n ) e iα n x-iβ ( j) n y
in the reflected and transmitted fields

∑ n∈Z (B j n , D j n ) e iα n x+iβ ( j) n y resp. ∑ n∈Z (A j n , C j n ) e iα n x-iβ ( j+1) n y) .
This is a linear operation between infinite vectors of the Rayleigh coefficients written as

{(B j n , D j n )} n∈Z = r j {(A j n ,C j n )} n∈Z , {(A j n , C j n )} n∈Z = t j {(A j n ,C j n )} n∈Z
with the linear reflection and transmission operators r j and t j , respectively. Similarly, the reflection and transmission operators r ′ j , t ′ j for illumination from below map the coefficient vector

{(B j n , D j n )} of the α-quasiperiodic incoming field ∑ n∈Z (B j n , D j n ) e iα n x+iβ ( j+1) n y
to the coefficient vectors of the reflected and transmitted fields

∑ n∈Z (A j n , C j n ) e iα n x-iβ ( j+1) n y) resp. ∑ n∈Z (B j n , D j n ) e iα n x+iβ ( j) n y . i.e. {(A j n , C j n )} n∈Z = r ′ j {(B j n , D j n )} n∈Z and {(B j n , D j n )} n∈Z = t ′ j {(B j n , D j n )} n∈Z .
Further, we assign to each profile Σ j an y-coordinate y j , for example y j = Y j (0) for a given parametrisation (X j (t),Y j (t)) of the profile Σ j , and define a diagonal operator γ γ γ j which maps a vector of pairs {(a n , b n )} n∈Z to the vector

{(a n , b j ) e iβ ( j) n (y j-1 -y j ) } n∈Z = γ γ γ j {(a n , b n )} n∈Z
If we introduce the infinite vector A A A 0 of the coefficients of the input wave

A A A 0 = {δ n0 (p z , q z )} n∈Z ,
then R 0 and T 0 are derived by the following marching procedure:

Set R M-1 = r M-1 , T M-1 = t M-1 ; Compute for j = M -1, . . . , 1 R j-1 = r j-1 + t ′ j-1 γ γ γ j R j ( I -γ γ γ j r ′ j-1 γ γ γ j R j ) -1 γ γ γ j t j-1 ; T j-1 = T j ( I -γ γ γ j r ′ j-1 γ γ γ j R j ) -1 γ γ γ j t j-1 ; Determine finally {(E + n , B + n )} n∈Z = R 0 A A A 0 , {(E - n , B - n )} n∈Z = T 0 A A A 0 .

Determination of the scattering matrices

For the application of the marching algorithm, one has to find finite-dimensional approximations of the scattering operators, i.e., scattering matrices, again denoted by r j , t j and r ′ j , t ′ j , for given j = 0, . . . , M -1. This means one-profile grating problems must be solved with incident waves from above and below for the profile Σ j . More precisely, one has to find the Rayleigh coefficients of the diffracted fields for input waves with z-components

( u + δ v + δ ) = ( 1 -δ δ ) e iα n x-iβ ( j) n y , ( u - δ v - δ ) = ( 1 -δ δ ) e iα n x+iβ ( j+1) n y , δ = 0, 1 . (12.107)
The choice of the indices n will be described in Sec. 12.6.1. First, we consider the calculation of the scattering matrices for a continuous interface Σ j . It separates two layers and the one-profile problem corresponds to the situation depicted in Figure 12.2. We denote the semi-infinite domains above and below the profile Σ = {(x, yy j ) : (x, y) ∈ Σ j } by G ± and by ε ± , µ ± the material coefficients above and below Σ, respectively. Thus, we keep the notation of Sec. 12.2.2, but the difference to the problem there is the occurrence of different incident waves from above and below and the fixed values ε 0 and µ 0 in condition (12.105).

For illumination from above, one has to solve the following problem: Setting

E z = { u + + u + δ in G + , u - in G -, B z (x, y) = { v + + v + δ in G + , v - in G -,
find α-quasi-periodic solutions of the Helmholtz equations

in G + ∆u + + ω 2 κ 2 + u + = ∆v + + ω 2 κ 2 + v + = 0 , (12.108) in G - ∆u -+ ω 2 κ 2 -u -= ∆v -+ ω 2 κ 2 -v -= 0 , ( 12.109)
where now κ 2 ± = ε ± µ ± -ε 0 µ 0 sin 2 ϕ . From equation (12.105) one gets the jump conditions on Σ

u -= u + + u + δ , ε -∂ n u - ε v κ 2 - - ε + ∂ n (u + + u + δ ) ε v κ 2 + = √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∂ t v -, v -= v + + v + δ , µ -∂ n v - µ v κ 2 - - µ + ∂ n (v + + v + δ ) µ v κ 2 + = - √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∂ t u -.
(12.110)

For illumination from below, we set

E z = { u + in G + , u -+ u - δ in G -, B z = { v + in G + , v -+ v - δ in G -.
The α-quasi-periodic functions u ± , v ± have to satisfy the Helmholtz equations (12.108), (12.109) and the transmission conditions

u -+ u - δ = u + , ε -∂ n (u -+ u - δ ) ε v κ 2 - - ε + ∂ n u + ε v κ 2 + = √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∂ t v + , v -+ v - δ = v + , µ -∂ n (v -+ v - δ ) µ v κ 2 - - µ + ∂ n v + µ v κ 2 + = - √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 + - 1 κ 2 - ) ∂ t u + .
(12.111)

Choosing as before u -, v -as single layer potentials (12.53), we derive from equations (12.110) and (12.111) the system of singular integral equations

( ε -κ 2 + ε + κ 2 - V + (I -L -) + (I + K + )V - ) w - √ ε v µ v √ ε 0 µ 0 ε + sin ϕ ( 1 - κ 2 + κ 2 - ) H + V -τ = U , ( µ -κ 2 + µ + κ 2 - V + (I -L -) + (I + K + )V - ) τ + √ µ v ε v √ ε 0 µ 0 µ + sin ϕ ( 1 - κ 2 + κ 2 - ) H + V -w = V , ( 12.112)
where the singular integral H + is defined by (12.30) with the fundamental solution Ψ ωκ + ,α . For illumination from above, the right-hand side is given by

U = -2 u + δ , V = -2 v + δ ,
whereas in the case of illumination from below

U = ε -κ 2 + ε + κ 2 - V + ∂ n u - δ -(I + K + )u - δ + √ ε v µ v √ ε 0 µ 0 ε + sin ϕ ( 1 - κ 2 + κ 2 - ) H + v - δ , V = µ -κ 2 + µ + κ 2 - V + ∂ n v - δ -(I + K + )v - δ - √ µ v ε v √ ε 0 µ 0 µ + sin ϕ ( 1 - κ 2 + κ 2 - ) H + u - δ .
In the case of a rod grating with a discontinuous profile, the domain G -is bounded. Using the single layer potential ansatz in G -, illumination from above is treated as before. Illumination from below can be treated by setting

E z = { u + + u - δ in G + , u - in G -, B z = { v + + v - δ in G + , v - in G -,
which results in the system (12.112) with the right-hand side

U = -2 u - δ , V = -2 v - δ .
Thus, in all considered cases the system (12.112) can be used to determine the scattering matrices. Moreover, it can be shown that the solvability of system (12.112) does not depend on ε 0 and µ 0 . Similar to the system (12.56), the equations are solvable if the ratios ε -/ε + and µ -/µ + do not belong to an interval on the negative axis. Thus, the applicability of the algorithm is independent of the incidence angles θ and ϕ as well as of the polarization.

Gratings with penetrating boundaries

In the following, we suppose that the interfaces Σ j are given by piecewise C 2 parametrizations

σ j (t) = (X j (t),Y j (t)), X j (t + 1) = X j (t) + d, Y j (t + 1) = Y j (t) , t ∈ R , ( 12.113)
i.e., the functions X j ,Y j are piecewise C 2 with

|σ ′ j (t)| = √ (X ′ j (t)) 2 + (Y ′ j (t)) 2 > 0 .
Moreover, the interfaces do not intersect, i.e. σ j (t 1 ) = σ k (t 2 ) only if j = k and t 1 -t 2 = dn. Additionally, we suppose that if a curve Σ j has corners, then the angles between adjacent tangents at the corners are strictly between 0 and 2π. To derive an integral formulation we rewrite the conical diffraction problem (12.8), (12.49), (12.11) using the notation

G M-1 G 2 G 1 G 0 G M E i y x d . . . . . n Σ M-1 Σ M-2 Σ 2 1 Σ 0 Σ
E z (x,y) = { u 0 + E i z u j , B z (x, y) = { v 0 + B i z in G 0 , v j in G j , j = 1, . . . , G M ,
with E i z = p z e i(αx-β y) , B i z = q z e i(αx-β y) . We seek α-quasiperiodic functions {u j , v j } N j=0 such that in G j ∆u j + ω 2 κ 2 j u j = ∆v j + ω 2 κ 2 j v j = 0 , (12.114) subject to the transmission conditions

on Σ 0          u 1 = u 0 + E i z , ε 1 ∂ n u 1 ε v κ 2 1 - ε 0 ∂ n (u 0 + E i z ) ε v κ 2 0 = √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 0 - 1 κ 2 1 ) ∂ t v 1 , v 1 = v 0 + B i z , µ 1 ∂ n v 1 µ v κ 2 1 - µ 0 ∂ n (v 0 + B i z ) µ v κ 2 0 = - √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 0 - 1 κ 2 1 ) ∂ t u 1 , (12.115)
and, for j = 1, . . . , M -1,

on Σ j          u j+1 = u j , ε j+1 ∂ n u j+1 ε v κ 2 j+1 - ε j ∂ n u j ε v κ 2 j = √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 j - 1 κ 2 j+1 ) ∂ t v j+1 , v j+1 = v j , µ j+1 ∂ n v j+1 µ v κ 2 j+1 - µ j ∂ n v j µ v κ 2 j = - √ ε 0 µ 0 ε v µ v sin ϕ ( 1 κ 2 j - 1 κ 2 j+1 )
∂ t u j+1 , (12.116) which satisfy the outgoing wave condition

(u 0 , v 0 )(x, y) = ∞ ∑ n=-∞ (E (0) n , B (0) n ) e i(α n x+β (0) n y)
for y > max

(x,t)∈Σ 0 t , (u M , v M )(x, y) = ∞ ∑ n=-∞ (E (M) n , B (M) n ) e i(α n x-β (M) n y)
for y < min (x,t)∈Σ M t .

(12.117)

The single and double layer potentials on one period Γ j = {σ j (t) : t ∈ [t 0 ,t 0 + 1]} of the interface Σ j corresponding to κ m are denoted by

S Γ j ,m φ(P) = 2 ∫ Γ j Ψ m,α (P -Q) φ(Q) dσ Q , D Γ j ,m φ(P) = 2 ∫ Γ j φ(Q)∂ n(Q) Ψ m,α (P -Q) dσ Q , with the α-quasiperiodic fundamental solution Ψ m,α = Ψ ωκ m ,α .
We present a recursive algorithm for solving (12.114 -12.117), which in each step treats a problem for one of the interfaces and therefore allows us to solve conical diffraction problems for gratings with an arbitrary number of layers on standard PCs (for the derivation, see App. B). The algorithm extends a method for in-plane diffraction, i.e., ϕ = 0, which was proposed by D. Maystre in Ref. 12.32 and described in detail in Ref. 12.4.

The starting point is to seek the solutions {u j , v j } M j=0 of (12.114-12.117) in the form

u 0 = 1 2 ( S Γ 0 ,0 ∂ n u 0 -D Γ 0 ,0 u 0 ) , v 0 = 1 2 ( S Γ 0 ,0 ∂ n v 0 -D Γ 0 ,0 v 0 ) , in G 0 , ( 12
.118)

u j = 1 2 ( S Γ j , j ∂ n u j -D Γ j , j u j ) + S Γ j-1 , j φ j-1 , v j = 1 2 ( S Γ j , j ∂ n v j -D Γ j , j v j ) + S Γ j-1 , j ψ j-1 ,      in G j , j = 1, . . . , M -1 (12.119) u M = S Γ M-1 ,M φ M-1 , v M = S Γ M-1 ,M ψ M-1 , in G M , ( 12.120)
with certain densities φ j , ψ j , j = 0, . . . , M -1. Again, the Helmholtz equations (12.50) and the outgoing wave condition (12.52) are satisfied. Note that the representations (12.118 -12.120) are unique, provided that the single layer potential operators V ( j) j-1 j-1 are invertible for j = 1, . . . , M, which will be assumed throughout.

The algorithm determines recursive relations

( φ j ψ j ) = Q j-1 ( φ j-1 ψ j-1
) , j = 1, . . . , M -1, (12.121) such that the functions {u j , v j } M j=0 fulfill the remaining transmission conditions (12.115) and (12.116). The initial densities (φ 0 , ψ 0 ) and the 2 × 2 operator matrices {Q j-1 } are obtained by the following scheme:

Introduce the coefficients

a j = ε j+1 κ 2 j ε j κ 2 j+1 , b j = µ j+1 κ 2 j µ j κ 2 j+1 , c j = √ ε v µ v √ ε 0 µ 0 ε + sin ϕ ( 1 - κ 2 j κ 2 j+1 ) , d j = √ µ v ε v √ ε 0 µ 0 µ + sin ϕ ( 1 - κ 2 j κ 2 j+1
) , (12.122) and determine Q j-1 by a backward recurrence for j = M -1, . . . , 1 as a solution of the operator equation ((

I + K ( j) j j -c j H ( j) j j d j H ( j) j j I + K ( j) j j ) A j - ( a j V ( j) j j 0 0 b j V ( j) j j ) B j ) Q j-1 = 2 ( V ( j) j j-1 0 0 V ( j) j j-1
)

. (12.123)
The initial values are

A M-1 = ( V (M) M-1M-1 0 0 V (M) M-1M-1
)

, B M-1 = ( L (M) M-1M-1 -I 0 0 L (M) M-1M-1 -I )
, (12.124) and the subsequent terms in (12.123) are derived by

A j-1 = ( V ( j) j-1 j-1 0 0 V ( j) j-1 j-1 ) - 1 2
(( K ( j) j-1 j -c j H ( j) j-1 j d j H ( j) j-1 j K ( j) j-1 j ) A j - ( a j V ( j) j-1 j 0 0 b j V ( j) j-1 j ) B j ) Q j-1 , (12.125) B j-1 = ( V ( j) j-1 j-1 0 0 V ( j) j-1 j-1 ) -1 ( I + K ( j) j-1 j-1 0 0 I + K ( j) j-1 j-1 ) A j-1 -2 ( I 0 0 I
) .

(12.126)

Having found A 0 and B 0 , the initial value (φ 0 , ψ 0 ) of (12.121) is a solution of the linear equation ((

I + K (0) 00 -c 0 H (0) 00 d 0 H (0) 00 I + K (0) 00 ) A 0 - ( a 0 V (0) 00 0 0 b 0 V (0) 00 ) B 0 ) ( φ 0 ψ 0 ) = -2 ( E i z B i z ) . ( 12

.127)

Then the solution above the grating is given by the integrals

u 0 = - 1 2 ( a 0 S Γ 0 ,0 ( I -L (1) 00 ) φ 0 + D Γ 0 ,0 V
(1) 00 φ 0 + c 0 S Γ 0 ,0 J

(1)

00 ψ 0 ) , v 0 = - 1 2 ( b 0 S Γ 0 ,0 ( I -L (1) 00 ) ψ 0 + D Γ 0 ,0 V (1) 00 ψ 0 -d 0 S Γ 0 ,0 J (1) 00 φ 0
) .

If desired, the field below the grating is found from the integrals

u M = S Γ M-1 ,M φ M-1 , v M = S Γ M-1 ,M ψ M-1 .
with the densities φ M-1 , ψ M-1 determined using the recursive relations (12.121).

Generalization of energy balance for lossy multilayer gratings

Resonance and non-resonance anomalies, differing in their nature, can be effectively explored in high-and low-conductive gratings, such as: surface plasmon excitations, Brewster and Bragg conditions, Rayleigh orders, groove shape and waveguide features, etc. Because of the s and p modes in conical diffraction being coupled through the boundary conditions, the associated problems are more general, and gratings can act as perfect absorbers and local-or/and surfacefield enhancers at any incidence polarization state. Knowledge of the accurate value of the absorption for a grating is very important for testing the correctness and reliability of developed programs. The energy balance is one of the basic accuracy criteria based on a single computation and it is generalized here in the case of lossy multilayer gratings. In this Subsection we derive important formulas for direct calculus of the absorption of multi-boundary gratings working in general conical mounts. Diffraction efficiencies for the reflected and transmitted orders in conical diffraction can easily be found from the corresponding Raleigh coefficients or boundary values, see (12.67)-(12.70).

If the multi-boundary grating is perfectly conducting, then for respective refractive indices ν 2 j = ε j µ j , Im ν M = ∞, and if there is no energy absorption in the grating layers, Im ν j = v j = 0, j = 1, . . . , M -1, then the energy conservation law is expressed by the standard energy criterion (see Ch. 2) under unitary normalization conditions

R = 1,
where R is the reflected energy.

If the grating is lossless, Im ν j = 0, j = 0, . . . , M, then the energy conservation law is expressed by a similar energy criterion

R + T = 1,
where T is the transmitted energy.

If Im ν j > 0 for some j = 1, . . . , M -1 and Im ν M = 0, then energy is absorbed in grating layers. Thus, the usual law of the energy conservation that the sum of efficiencies of all reflected and transmitted orders should be equal to the power of the incident wave, does not hold. In a general case,

A + R + T = 1, (12.128)
where (see (12.157))

A = 1 β Im ∫ Γ 0 ( ε 0 ε v ∂ n E z E z + µ 0 µ v β ∂ n B z B z ) - κ 2 0 β κ 2 M Im ∫ Γ M-1 ( ε M ε v ∂ n E z E z + µ M µ v ∂ n B z B z
) is called the absorption coefficient or simply the absorption in the given multilayer diffraction problem. If also Im ν M > 0, then T = 0 and it holds

A + R = 1 (12.129)
with the absorption (12.156)

A = 1 β Im ∫ Γ 0 ( ε 0 ε v ∂ n E z E z + µ 0 µ v ∂ n B z B z
) .

The requirements (12.128), (12.129) are convenient single computation tools to check the quality of the numerical solution. Besides being physically meaningful, expression (12.157) is very useful as one of numerical accuracy tests for computational codes and especially important for x-ray-EUV gratings, photonic crystals, metamaterials, and perfect absorbers where absorption plays a predominant role. In the lossy multilayer case, one needs an independently calculated quantity A to verify (12.128). For such a quantity, we use the absorption integrals defined in Ref. 12.33 and derived in Appendix C using the second Green formula and integration by parts. The expressions derived from the boundary integral equation theory are important for calculating the absorption of general multi-boundary gratings working in any diffraction mount at any polarization state. The boundary absorption integrals developed and tested have been found to be an accurate and universal tool for calculating of the energy balance of various periodical structures having separated or penetrating boundaries. The results of absorption calculus of a bare metallic grating with shallow grooves, photonic crystal supporting polariton-plasmon excitation and x-ray-grazing-conical-diffraction multilayer grating are demonstrated in Sec. 12.9.

Remark 12.5.1 A generalization of the energy balance presented for a multilayer absorption grating in classical and conical diffraction is based on computations of the respective absorption integrals by values of the field and its derivatives on a boundary. It has not only intuitive significance but the same rigor, namely in the sense of generalized functions or distributions, and way to deduce as more simple energy criterions for perfectly conducting and lossless gratings (see in Ch. 2). A derivation of expressions considered for finding the absorption quantity A as well as the interpretation of the results obtained bear only on Maxwell's equations or Green's theorem and boundary conditions. The computation of A itself is not connected with a specific rigorous method which is used for near-zone field calculus. Thus, the present in-plane and off-plane energy balance generalizations for multilayer absorbing gratings can be considered as much universal and useful as well known energy conservation laws for perfectly conducting and lossless gratings.

Implementation and algorithmic enhancements of multilayer solvers

To handle effectively various grating types, the different multilayer schemes can be used to solve respective diffraction problems, i.e. Penetrating or Separating solvers. The Penetrating solver described above is more general, since it allows the y-projections of the boundaries to be overlapping that is vital in the modeling of many coated gratings. However, when the grating boundary profiles are strictly separated, the problem (12.104)-(12.106) can be treated using certain robust algorithms for modeling layered gratings. Therefore, the Separating solver based on the S-matrix multilayer algorithm can be, for particular problems, several times faster and more accurate than the Penetrating one.

There are three basic sources of numerical errors arising from an integral equation implementation: (i) replacement of an integral equation by a finite system of linear algebraic equations; (ii) inexact evaluation of matrix elements; (iii) inaccuracy of solution of the linear system. For errors of type i, in many cases a priori estimations via functional-analytic and operator-theoretical methods are available, which, at least, can moderate one's optimistic expectations about the overall convergence rate. Combinations of moments, Galerkin, collocation, and fully discrete methods with balanced convergence properties are known as numerical schemes of discretization [12.22]. Errors of type ii in these methods are commonly attributed to numerical quadratures. In periodic diffraction problems, in contrast to diffraction on a compact obstacle, there is one more source of ii-type numerical errors: evaluation of lattice Green functions and their derivatives. The problem is seen from the well-known kernel functions representation described above. This is the most difficult error type arising from solving of grating-like diffraction problems and particularly at small λ /d ratios. Errors of type iii, as well as direct round-off errors, are negligible in most cases provided the numerical scheme in use is stable and the problem "generic". That is also true for iterative linear system solvers used in our codes, like GMRES-or FOM-based software (see Sec. 12.4.3).

In order to reduce time for computation matrices of the above operator equations, two further basic enhancements (cache for exponential functions (plane waves described in Sec. 12.4.6) at the algorithmic level are used in our codes: cache for kernel functions, and cash for repeating pairs or quads of layers. They assume a large time-memory trade-off. The amount of memory required for cache can be calculated in advance in each case and adjustments (cache off or partial) are done automatically. More acceleration can be reached in some cases, e.g. if one uses iterative algorithms to solve a linear system of algebraic equations (see Sec. 12.4.3).

Implementation of multilayer schemes

Here we describe an implementation of the S-matrix algorithm combined effectively with the conical integral equations formulated for solving such multilayer grating problems. We discuss briefly the numerical solution of the system (12.112). As mentioned before, the scattering matrices are obtained by solving one-profile equations with a finite number of illuminations (12.107). The indices n of these illuminations should be chosen such that additionally to the diffracted outgoing modes the Rayleigh coefficients of some evanescent modes are also taken into account. Let the grating formed by the profiles Σ j , which separate optical materials with the parameters ε j , µ j and ε j+1 , µ j+1 . The indices of propagating modes are characterized by the values n such that β ( j) n ≥ 0 above Σ j and β ( j+1) n ≥ 0 below Σ j . Suppose that their number is P u ≥ 0 above and P d ≥ 0 below the profile. Further, we fix numbers k u and k d of evanescent modes which are important to keep in the scattering matrices. This results in quadratic reflection matrices r j and r ′ j of order 2(P u + k u ) × 2(P u + k u ) and 2(P d + k d ) × 2(P d + k d ) for illumination from above and below, respectively. The transmission matrices t j and t ′ j are rectangular of dimension 2(P d + k d ) × 2(P u + k u ) and 2(P u + k u ) × 2(P d + k d ) for illumination from above and below, resp.

These matrices are constructed columnwise from the scattering amplitudes of the solutions of the equation (12.112) with right-hand sides of index n within the fixed range. Note, one has only once to discretize the integral operators in (12.112) and the LU-decomposition of this discrete matrix provides the solutions immediately and, hence, all four scattering matrices simultaneously. Therefore, we use a direct solver with LU-decomposition for computing the scattering matrices. It should be noted that modern implementations of the LAPACK and BLAS software packages on multi-processor/core/thread machines makes direct solving a competitive alternative to iterative solution methods even for very large systems, N 10000. Expressions (12.149) and (12.150) allow us to find amplitude matrices by a recursive procedure beginning with the lower medium. To do this, we have to know, in a general case, four matrices of scattering amplitudes and perform two matrix inversions in each iteration step. The computation time for one-boundary problems was shown to scale quadratically with the main accuracy parameter N (see Fig. 12.5). The computation time is also linearly proportional to the number of boundaries. The memory cache for amplitude matrices of multi-layer grating problems (e.g. photonic crystals) with the same boundary profiles and the same pairs or quads of layers can be used (see Section below).

Efficient implementation of the penetrating solver should use the modern implementations of the LAPACK and BLAS software packages and their analogues on multi-processor/core/thread machines. Although the algorithm requires a larger number of matrix-matrix multiplications compared to the separating solver and even the inversion of dicretization matrices, even quite complicated problems can be solved on a modern PC in reasonable time.

Cache for kernel functions

Matrix elements of discretized integral equations are kernel functions of one of four types considered: single-layer potentials, double-layer potentials, normal derivatives of single-layer potentials, and tangential derivatives of single-layer potentials. Any of these kernel functions for the given layer has two vector arguments: the source position x and the observation point x 0 . The value of the kernel function depends on the difference vector d = xx 0 .

There are a number of cases of practical interest when the same difference vector d corresponds to more than one pair {x, x 0 }. Typical situations include:

• conformal layers; upper and lower boundaries of such a layer are obtained from each other by a vertical shift;

• more generally, layers whose boundaries are congruent by a translation (not necessarily in the vertical direction);

• rectilinear segments of boundaries, if collocation points are uniformly distributed along such a segment.

In all these situations, it is possible to re-use values of kernel functions calculated earlier. The program stores the data: type of potential and difference of arguments vector (∆x, ∆y) -in lexicographical order. Fast search and insertion are provided by a binary tree structure [12.34]. Memory expenditure for the kernel function cache is of the order cN 2 per layer, where N is the maximum number of discretization points on the boundaries, and the constant c incorporates the size of data structure corresponding to each node of the tree. If no further layer has a refractive index of the current layer, then the cache gets overwritten as the solver proceeds to a new layer. However, it is quite typical to have a multilayer structure with repeating indexes, in which case the kernel function computed for one layer has a chance to be re-used in another layer. Note that the constant c is less the more effective the cache is (that is, the more repetitions occur). To save memory, single precision values are used for the difference components ∆x, ∆y. This approach does not compromise accuracy to any noticeable extent.

Cache for repeating pairs or quads of layers

The memory cache for scattering amplitude matrices (computation matrices of the operator equations considered) of multilayer grating problems with separating boundaries with the same boundary profiles and the same pairs or quads of layers can be used. The actual number of identical pairs or quads of layers can be large, up to a thousand for hard x-ray grating applications. For flexibility and possible reuse of scattering matrices of the Separating solver in multi-stack grating structures with repeated layer patterns, the dynamic caching procedure using a hash function for fast storing and extracting of boundary and adjusting layer basic parameters is initialized separately for each interface starting from the bottom. In such a procedure previously calculated instances are taken for reuse in accordance with hash function values. Computation time of the efficiency of a coated grating with many of repeated pairs or quads of layers and equal separated boundaries can be decreased by orders of magnitude using such memory cache for scattering matrices. For example, if the number of pairs of layer with the same boundary is more than 50, then the computation time for typical x-ray grating problems decreases more than 100 times-see Fig. 12.6.

12.7 Modifications of integral methods for very small wavelength-to-period ratios It is well-known that solution of the 2D Helmholtz equation with any rigorous numerical code meets with difficulties at small λ /d. While the standard IMs are robust, reliable and efficient, they exhibit poor convergence and loss of accuracy in the high-frequency range due to numerical contamination in quadratures. Increasing matrix size and enhancing computation precision, as well as application of convergence speed-up techniques, which are successfully explored in low-and mid-frequency ranges, lead to unreasonably stringent requirements for computing times and storage capacities in high and, especially, ultrahigh frequency ranges (d/λ > 10 3 ). For various kinds of integral equations and approximation technique used for solving diffraction grating problems the computation accuracy is mostly determined from the accuracy of computing the fundamental solution (see, e.g., Sec. 12.4). In order to approximate a wave with wavenumbers k j = ν j d/λ in accordance with the Rayleigh criterion, in the Helmholtz equation one needs to use about 10 (usually from 5 to 20; it depends on a groove profile) discretization points per wavelength. So, for very large wavenumbers, say k j ∼ 1000, discretization matrix size should be of the order N ∼ 10 4 , a huge enough number even for modern PCs. The inaccuracy of computation of the fundamental solution, together with some rounding errors, increases significantly, up to totally divergent results, if one goes far from this simple rule of thumb. The term modified integral method used in publications with regard to the PCGrate software was introduced with flexibility in mind. More precisely, it is meant to be '"modifiable" or "tunable", however we keep the earlier term as a label. In a narrow technical sense, the MIM is characterized in this Section by a number of simple modifications required for the standard IM, similar to the one described in Ch. 4, to transform it to the MIM, together with relevant discussions.

The boundary integral equation theory is so flexible and complex that we can point out a few areas of its modifiability. (1) In the physical model, one can choose boundary types (periodical or non-periodical, closed or non-closed, smooth or having edges, randomly rough or deterministic, etc) and boundary conditions (rigorous or non-rigorous, perfect or finite conducting, extending, etc). (2) In the mathematical structure, integral representations using various potential operators and/or integral formulas together with multi-layer schemes can be considered. (3) In the method of approximation and discretization, numerical scheme of discretization (Nyström, or collocation, or Galerkin, or Method of Moments, hybrid, etc), basis and test functions (piecewise constant, or trigonometric polynomials, or splines, or delta, or Lagrange polynomials, etc), and including treatment of corners in boundary profile curves can be chosen. [START_REF] Goray | Spectral separation of the efficiencies of the inside and outside orders of soft-x-ray-extreme-ultraviolet gratings at near normal incidence[END_REF] In the low-level details, one can define methods of calculations of kernels (direct methods using Hankel or exponential functions, or Ewald's method, or high-order summations, etc); and using acceleration techniques (Kummer or Euler-Knopp summation, or single-term corrections, etc), meshes of sampling (collocation) points (uniform or non-uniform), quadrature rules (trapezium, or Gaussian, or more sophisticated); solutions of linear systems (direct methods or iterative solvers). [START_REF] Goray | Application of the boundary integral equation method to very small wavelength-to-period diffraction problems[END_REF] In the implementation improvements, one can use caching of repeating quantities (exponential and kernel functions, pairs or quads of layers in multilayer structures, etc), treatment of Rayleigh orders, etc. In this Section, special attention is paid to important aspects of the presented MIM for small λ /d diffraction problems in connection with (3)-( 5), and also, briefly, to some details of the numerical implementation. More about the MIM im-plementation can be found in the documentation devoted to the PCGrate software and also in specific references [12.1].

Diffraction from 1D multilayer gratings with arbitrary boundary profile shapes, including edges, is considered in this Section in a general case of off-plane mounts. The term 1D multilayer refers to a general grating on a planar surface of arbitrary conductivity which is periodic in one direction, constant in the second, and has a finite number of borders and layers in the third. The actual number of identical or different borders and layers can be large enough, up to a few thousand for hard x-ray grating applications. Though various approximated analyses developed for the treatment of such challenging diffraction problems enjoy more or less successful application [12.35], they are always plagued with uncertainties which make comparison between rigorous and non-rigorous approaches difficult. In the present study, special attention is paid to all aspects of the MIM for λ /d ≪ 1 ratios. A few commercial and non-commercial solvers based on the MIM are used in this Section and also in Sec. 12.9 to analyze the diffractive properties of various bulk and multilayer gratings including those having real boundary profiles of the polygonal type obtained by averaging measured data from Atomic Force Microscopy (AFM).

Approximations

As to the multilayer schemes implemented, there are no substantial differences between the well established approaches suitable for resonance domains (see Secs. 12.5 and 12.6) and the MIM in these higher levels of the multilayer boundary integral equation theory. We use both the Penetrating and Separating solvers to treat efficiencies of multilayer x-ray-EUV gratings having many boundaries with thin structure including random micro-and nano-roughness (see Sec. 12.8 and Sec. 12.9.9). However, the mid-and low-level MIM structures including the method of discretization have a few important peculiarities described below.

It is well-known the convergence and accuracy of IMs depend greatly on an appropriate choice of the discretization scheme and respective quadrature method for solving integral equation systems. As a rule, a Nyström discretization or a collocation method, as well as a Method of Moments or a Raleigh-Ritz-Galerkin approach, which are not described here, or their combination, is a good choice to treat both general and particular diffraction problems. The sampling points of unknown functions can be distributed on some uniform or multi-scale grids. In lowand mid-frequency ranges, better results are often obtained using equidistant steps along the arc length. Another possible function of the distance between collocation points is prescribed by equal steps along the x-axis perpendicular to the grooves.

In the MIM, the fastest Nyström method with the rectangular quadrature rule is used (see Sec. 12.4). Such a simple, fully discrete method combined with some matrix element modifications works well for shallow smooth boundary profiles and, particularly, at small λ /d. In the presence of a boundary profile with corners (piecewise linear), another approach can be effective. The sampling and quadrature nodes are set in such a way that every corner lies halfway between the adjacent nodes and no curvature-like single-term corrections (see (12.82)) are added to diagonal matrix elements. Let us match a solution in N midpoints σ

(t i+1/2 ) of [σ (t i ), σ (t i+1 )] by setting φ -(σ (t i+1/2 )) = (φ -(σ (t i+1 )) + φ -(σ (t i )))/2.
Then we obtain a linear system of algebraic equations for φ -(σ (t k )), k = 1, . . . , N similar to (12.80)

φ -(σ (t i+1/2 )) + N ∑ k=1 c i+1/2,k φ -(σ (t k )) = b(σ (t i+1/2 )), σ (t k ) = (X(t k ),Y (t k )) ∈ Γ . (12.130)
In this approach the period of integration is divided by a number of segments equal to the number of corners on the boundary profile. Thus, the sampling points and the quadrature nodes can be put at same locations, as in (12.80), or interlacing by a half-segment shift, as in (12.130). The choice between two these complementary approaches depends on desired accuracy of computations and time requirements. For shallow boundaries with a thin structure (multi-polygonal) including roughnesses, the approach of (12.130) may have faster convergence. However, for boundaries with several rather long and abrupt slopes the approach of (12.80) may be preferable, especially if one uses (1) the curvature single-term correction by adding the corner term to diagonal matrix elements or (2) the mesh grading together with the appropriate quadrature rule, as in a case of deep gratings.

Instead of the direct summation algorithm used in the MIM and also in the IM of Ch. 4, more sophisticated methods can be implemented to accelerate the computation of the integral equation kernels, like as Ewald's methods (see Sec. 12.4.5.2). Unfortunately, it has turned out numerically that such approaches, at least those known for us, are not efficient for very small λ /d diffraction grating problems. Thus, the MIM in a narrow sense is an approximation approach with a simple discretization that also specifies a summation method for kernels.

Convergence and accuracy with and without speed-up technique

It is well known that the number of discretization points per wavelength used in the various IMs can be reduced significantly, up to an order of magnitude, when λ /d and H/d become small. The question is how small it might be for very small λ /d diffraction problems, say for λ /d 10 -3 . To accelerate convergence of the series representing kernels, different acceleration techniques can be applied (see Sec. 12.4 and also Ch. 4). In Figs. 12.7-12.10, convergence of the IM is demonstrated for an analytical case of diffraction from a plane transmission interface prescribed by a zero-depth sinusoidal profile at normal incidence in a vacuum with the lowermedium refractive index of ν 1 = 1.5 and for different λ /d. Note that for all numerical examples in this Subsection, the number of positive and negative terms accounting in kernel functions was chosen at N/2 (see in Section 12.7.3).

For λ /d = 1 in Fig. 12.7, the convergence rate reached with speed-up techniques (all single-term corrections in kernels are used) is high, with the energy balance error of ∼ 10 -6 in both polarization states for the number of discretization points N = 10. For λ /d = 10 -1 in Fig. 12.8, the convergence rate reached with speed-up techniques is medium, with the energy balance error of ∼ 10 -5 in both polarizations for N = 100. For λ /d = 10 -2 in Fig. 12.9, the convergence rate, again obtained with speed-up techniques, is low, with the energy balance and transmitted energy errors of ∼ 10 -3 in both polarizations for N = 500. The difference between the TE and TM transmitted energies for N < 300 is seen to be large, ∼ 10 -1 . For λ /d = 10 -3 in Fig. 12.10, the convergence rate calculated with speed-up techniques is very low, with the Energy balance error of ∼ 10 -2 in both polarizations for N = 10 3 . As seen from the figure, the convergence of the series deteriorates for N > 1000 as the distance between the kernel function's arguments tends to zero. In contrast to the plots of Fig. 12.10, the results for λ /d = 10 -6 obtained without application of any speed-up techniques exhibit the fastest convergence rate, with a negligible energy balance error of ∼ 10 -16 for N = 2 only, and are equivalent to analytical calculations. In Fig. 12.11, convergence of the standard integral method is demonstrated in respect to the main cut-off parameter N for another analytically amenable case, i.e. of x-ray diffraction from a plane absorbing interface (grazing incidence in vacuum of non-polarized radiation to plane Au surface prescribed by a zero-depth sinusoidal profile) for λ = 1 nm, θ = 89 • , and different λ /d. The refractive indices of Au for all examples in this Section were taken from Ref. 12.36. For λ /d = 10 -2 , the convergence rate reached using speed-up techniques (i.e. by the standard IM) is high, with the reflected energy error of ∼ 4.9 × 10 -6 for the number of discretization points N = 40 (the exact reflectance value is 0.7999). For λ /d = 10 -3 , the convergence rate reached with speed-up techniques is medium, with the reflected energy error of ∼ 10 -3 for N = 200. For λ /d = 10 -4 , the convergence rate, again obtained with speed-up techniques, is low, with the reflected energy error of ∼ 6.2 × 10 -3 for N = 10 3 . For λ /d = 10 -5 , the convergence rate calculated with speed-up techniques is very low, with the reflected energy error of ∼ 7.7 × 10 -2 for N = 2 × 10 3 . In contrast to the plots of Fig. 12.11, the results for extremely low λ /d of 10 -7 obtained by the MIM without application of any speed-up techniques exhibit the fastest convergence rate with a negligible reflected energy error of ∼ 10 -16 for N = 2 only and are equivalent to analytical calculations. Thus, we see for this grazingincidence absorbing example the same behavior of kernel functions as in the previous absolutely different case of the normal incidence on the lossless medium.

As one can see from Figs. 12.7-12.11, at least one discretization point per wavelength is required to reach efficiency convergence for the standard IM. In contrast to that, the MIM with the simple, however very important, changes in respect to the described above standard IM, i.e., without applying acceleration convergence terms at low λ /d only, works accurately and ultra-rapidly despite the very small number of discretization points per wavelength used in the approach. For example, if a period includes N = 10 2 and λ /d = 10 -4 , there is only 10 -2 point per wavelength required for the MIM. While the results presented in Figs. 12.7-12.11 may certainly be different for various realizations of the IMs and of the speed-up techniques used, the overall pattern remains the same. Shallow gratings and rough mirrors exhibit a similar behavior for very small λ or λ /d (very large wavenumbers k) in the x-ray and EUV ranges (see Figs. 12.12 and 12.13). In this case, however, the effective boundary profile depth ∼ H cos θ , the bilayer thickness ∼ λ /(2 cos θ ′ ) (θ ′ = arcsin[(sin 2 θ cos 2 ϕ + sin 2 ϕ ) 0.5 ]), and the effective radiation wavelength ∼ λ /(ν j cos ϕ j ) must be of the same order of magnitude. In the present approach, the peculiarity described in Ref. [START_REF] Bremer | On the Nystrom discretization of integral equations on planar curves with corners[END_REF] ("Introducing known speed-up terms in integral methods produces an adverse numerical effect because of the ensuing uncontrolled growth of coefficients in analytically (or numerically-Goray & Schmidt) improved asymptotic estimations") takes into account mostly for the case of shallow x-ray-EUV gratings working at very small λ /d and including, if any, random roughness (for more, see Remark below and also Section 12.8). As shown in this Subsection with all speed-up options turned off, it is possible to obtain for the most difficult problems of small λ /d ratios surprisingly rapid convergence, and an energy balance very close to 1. The most important among the convergence speed-up options which have to be switched off in this case is the allowance for logarithmic singularity, and second in importance is the correction applied to account for the cut-off terms in the expansions of kernels (the Aitken's δ 2 single-term correction in our case (see Sec. 12.4.5.1)). Switching off the curvature single-term correction is of lower but not minor significance on the way to reaching fast convergence. Such calculations at very low λ /d also depend significantly on the actual summation rule chosen for the kernel functions that will be discussed in next Subsection.

Remark 12.7.1 The same rule for the relations between basic grating and light parameters and reaching the maximum diffraction efficiency in a desired order is, on the whole, valid for longer wavelengths, too. For example, the MIM with speed-up options turned off can be applied also for echelle gratings working at very high orders (very low λ /d) and H cos θ /d ≪ 1 [12.5, 12.35]. Thus, the record of rigorous computations was achieved for the r-10 EXES echelle grating of the NASA project SOFIA with d = 7.62 mm working in the -1431 order at a wavelength

Summation rules for kernel functions and energy balance

For many practical cases, there are no big problems to reach fast convergence and sufficient accuracy of results obtained by only varying the major accuracy parameter N j , which is usually the same for all boundaries of gratings layers: N j = N. The N values of 100-400 provide good accuracy commonly, with the exception of the following difficult cases: very deep (in respect to period and/or wavelength) boundaries; real boundary profiles with super fine structures including random roughness; very close boundaries; extremely grazing incidence; bad points on Rayleigh wavelengths, resonance anomalies of different kinds; high order echelles; high conductivity (especially for the TM polarization); some others and, especially, a combination of a few of these cases. For such complex problems, an increase in the number of discretization points may become necessary. However, to obtain accurate data for hard examples of computations, i.e. at very low λ /d, optimization of another accuracy parameter should be fulfilled.

In addition to N, there is one more important code parameter, namely the "Maximal number of accountable plus or minus terms" that describes a number of positive and negative terms accounting in kernel functions. This is the number of grating adjacent periods accountable in expansions of Green functions and their derivatives due to the quasi-periodicity property of the fields. In the simplest case typical of real problems, all kernels are truncated symmetrically in respect to the upper and lower regions and equally for any j-th boundary:

P± j = Pj = P ≈ gN j = gN. ( 12.131)
The "truncation ratio" g is optimized at small values of N and is kept constant as N increases. It has been found [12.25] that g = 1/2 is a reasonably good choice for most practical compu- tations, and especially in the short wavelength range. Typical dependencies on P for the above example with λ /d = 10 -2 are shown in Fig. 12.14. The energy balance is closer to 1 in both polarizations and TE/TM transmitted energies are close to each other at P = 50% of N, with divergence seen to set in at smaller and larger values of g. For another, very different, example of the absorbing x-ray grating working at grazing incidence, one can see similar dependencies of the energy balance on P in Figs. 12.15 and 12.16 for different numbers of discretisation points. The energy balance is close to 1 for both integral methods considered at P ≈ 50% of N = 10 2 with a very high rate of convergence for the MIM and very slow convergence for the standard IM, similar to the convergence dependencies on N presented above. While the MIM has the long-range of high accuracy converged results from P ≈ 40% to P ≈ 70% in Fig. 12.15, only two points in a curve for the standard IM have the energy balance values close to 1, with divergence seen at both sides from these points. Similar behavior is seen in Fig. 12.16 for the high value of N = 400, where again the energy balance is close to 1 for the MIM and the standard IM at g ≈ 0.5.

While today this rule is no more than empirical, there can be no doubt whatsoever that this choice is valid, and this has been verified in many realistic examples during the past two decades. Note that in the IM developed by D. Maystre during the later 70s [12.37], g = 2/3 for the resonance domain (λ ∼ d) and should be varied for different λ /d . It is worth noting that g = 2/3 is worse than g = 1/2 because the computation time is proportional to 2 PN 2 . It is interesting that the first "good" point in Fig. 12.15 for the standard IM is close to the value of g = 0.5, i.e. our "golden rule", and the second "good" point-to the value of g = 0.7, which agrees well with the rule of g = 2/3 given earlier for the standard IM. The present golden rule is also approximately satisfied for the all examples of numerical results given in Ch. 4.

To reduce computing time for matrices of the discretized operator equations, a few enhancements at the algorithmic level are used in the MIM: cache for kernel functions, cache for exponential functions, and cache for repeating pairs or quads of layers of multilayer gratings (see Secs. 12.4.6, 12.6.2, and 12.6.3). They assume a big time-memory trade-off at small λ /d. The amount of memory required for cache can be calculated in advance in each case and ad-justments (cache off or partial) are done automatically. The computation algorithm group in PCGrate codes enables one also to choose an algorithm for solving linear systems of algebraic equations. It can be either the direct Gauss or the non-direct FOM method (see Sec. 12.4.3). Note that in the Penetrating solver, linear systems are solved with the Gauss algorithm alone.

Remark 12.7.2 One more important note regarding the energy balance summation for very small λ /d problems appears to be pertinent here. The Green function and its derivative members tend to big values near Raleigh wavelengths when the y-component β (0),(M) n of the n-th diffraction order wave vectors in the upper medium or/and in the lower medium (for transmission gratings) tends to zero (see (12.20)). This means that the diffraction order becomes grazing or even close to evanescent. Its efficiency may be rather high from the physical point of view or/and diverge from the mathematical point of view (it depends also on N). It is well known from the diffraction theory that the efficiency of strictly grazing propagating, as well as of all evanescent, orders is zero. Moreover, various rigorous and approximate methods valid for shallow gratings operating at small λ /d, as well as all experimental data suggest convincingly that the efficiency decreases rapidly with increasing modulus of the diffraction order number. As a rule, the efficiencies of such grazing orders should be very close to zero and much less than the inaccuracy of computations. Thus, rather big and diverging efficiencies of high number grazing orders should be excluded from the energy balance consideration, for example, starting from a high order which becomes increasing in efficiency.

Analysis of rough gratings using quasi-periodicity and Monte Carlo calculus

Multi-wave and multiple diffraction, refraction, absorption, waveguiding and wave deformation govern, to a considerable extent, scattering of x-ray and EUV radiation and cold neutrons from nanoroughness of continuous media. Inclusion of these pure dynamic effects, which requires application of electromagnetic theory, permits one to calculate the absolute intensity of coherent (specular or diffraction order) components and describe adequately the intensity distribution of the non-coherent (diffuse) components which may have resonance peaks. Some surfaces are deterministic, e.g., perfect gratings, and some are random, e.g., polished mirrors). Some surfaces are 1D, e.g., one-periodic (classical) gratings and cutting mirrors, but most are 2D, e.g., bi-periodic gratings (bigratings), ocean surfaces, and surfaces with atomic scale roughness. Any number of possible combinations between these four characteristics may be present in real structures, e.g., 1D deterministic gratings modulated with 2D random roughness. Despite the impressive progress reached recently in development of exact numerical methods of investigating wave diffraction from boundary roughness [12.38, 12.39], the present authors are aware only of asymptotic and perturbation approaches to the analysis of x-ray and neutron scattering for 1D and 2D rough surfaces, such as the scalar Kirchhoff integral, parabolic wave equation, Rayleigh method, Born approximation, distorted-wave Born approximation, and a few others [12.40, 12.41]. The MIM and other rigorous approaches identified that the intensities of xray-EUV scattering at boundaries with random roughnesses may differ considerably (by a few times) from the values derived with the use of various approximate models [12.5, 12.6]. It was found that the MIM operates equally well with nano-roughness of any kind and shape which obey arbitrary statistics of distribution (not necessarily periodic or Gaussian, or fractal, etc.).

There are two classical and equivalent approaches, with some restrictions in each of them, to model rigorously randomly-rough 1D and 2D surfaces. The most general and timeconsuming one is to use large surface lengths of many wavelengths. In this approach some window functions and tapered (narrowing) beams can be used to restrict the illuminated range and avoid numerical difficulties at endpoints. The second widely-explored approach is to use periodic boundary conditions (quasi-periodicity of Floquet-Bloch modes). This method uses an infinite beam (plane wave) and assumes that the random rough surface lengths repeats itself for given large periods having some numbers of random asperities. That means using infinite grating samples together with intensive Monte-Carlo simulations. Examples of the both famous approaches are well described in the literature, see e.g. in Refs. 12.38, 12.42-12.44. From the theoretical and numerical reasons we thought it convenient to use the large-period-grating model to analyze shallow randomly-rough gratings in the x-ray-EUV range. This classical model for computation of bulk or few-border rough mirrors and quasi-gratings is applied in PCGrates and other of our codes to calculate multilayer rough mirrors and gratings, as well as multiple quantum dot or quantum molecular ensembles with most realistic border profiles having irregularities of any kinds, including real ones, i.e. measured by AFM, Transmission Electron Microscopy, Near-field Scanning Optical Microscopy, etc, or derived from simulations using a continuum growth model of multi-scale reliefs [12.45-12.47].

Diffraction from 1D surface grating-like structures with shallow boundary profile shapes is considered in this Section for the sake of simplicity for bulk gratings working in conical diffraction at small λ /d ratios. A generalization to a multilayer case is straightforward. The integral equations developed in the previous Sections are used in the present Section to analyze the diffractive properties of bulk gratings with real-profile boundaries having random roughnesses. The Section also reports on the electromagnetic solution of reflection from 2D rough surfaces in short waves using boundary integral equations for gratings in conical diffraction and Monte Carlo simulations. The general equivalence rule for determination of the efficiencies of reflected orders of bigratings (2D) from those calculated for classical (1D) gratings is derived.

Scattering intensity, absorption, and energy balance of rough 1D gratings

For a given incident plane wave with wave vector k = (α, -β , γ) = k + (sin θ cos ϕ ,cos θ cos ϕ , sin ϕ ), the reflected and transmitted diffraction orders of number n have the wave vectors

k ± n = (α n , ±β ± n , γ) = k ± (sin θ ± n cos ϕ ± , ± cos θ ± n cos ϕ ± , sin ϕ ± ), with (k ± ) 2 -γ 2 = α n 2 + (β ± n ) 2 , (β ± n ) 2 ≥ 0. Since the z-dependence of all functions is given by exp(iγz) tan θ ± n = α n /β ± n , ϕ + = -ϕ , ϕ -= arcsin(k + sin ϕ /k -).
By convention, the outgoing angles θ ± n of the reflected and transmitted orders (to ensure that θ + 0 = -θ ) are taken from the interval [-π/2 , π/2], as well as ϕ + and ϕ -. The p-and s-components of the E-fields of the plane waves (incident and diffracted) are defined with respect to the grating normal n = (0, 1, 0). We define the vectors s orthogonal to the plane spanned by k and the grating normal and p lying in that plane (see Sec. 12.2):

s = k × (0, 1, 0)/|k × (0, 1, 0)| , p = s × k/|k|.
If k = (0, k, 0), we set s = (0, 0, 1) and hence p = (1, 0, 0). Since for a plane wave the electric field E is orthogonal to the wave vector, (E, k) = 0, one can decompose E E = (E, s) s + (E, p) p.

The polarization angles of the wave are defined as

δ = arctan[|(E, s)|/|(E, p)|], ψ = -arg[(E, s)/(E, p)],
where δ ∈ [0, π/2], ψ ∈ (-π, π]. Such a representation of polarization angles is useful to define polarization states and polarization properties of incoming and diffracted waves in diffraction grating applications, see, e.g., Examples in Sec. 12.9.

The efficiency of a diffracted order represents the proportion of power radiated in each order. Defining the power as the flux of the Pointing vector modulus |P| = Re (E × H)/2 through a normalized rectangle parallel to the (x, z)-plane, the ratio of the power of a reflected or transmitted propagating order and of the incident wave gives the conical diffraction efficiency η ± n of this order in the simple form (see (12.67), (12.68)). For the reflected orders we have

η + n = β + n β ( ε + ε v |E + n | 2 + µ + µ v |B + n | 2 ) ,
where the formulas for E ± n , B ± n are given by (12.69) and (12.70). If Im k -> 0 then there are no transmitted orders. Thus, the usual law of energy conservation, that the sum of efficiencies of all reflected and transmitted orders should be equal to the power of the incident wave, does not hold. Instead, some part of the power is absorbed in the substrate. If the grating is absorbing, then conservation of energy is expressed by a criterion

R + A = ∑ β + n >0 η + n + A = 1, (12.132)
where R is the sum of the reflection order efficiencies and A is the absorption in the singleboundary off-plane problem that can be computed from integrals of the solution of the partial differential formulation of conical diffraction (see (12.74)). For the general elliptically polarized incident light in conical diffraction, the reflected efficiency can be found as

η + n = |C + n (θ , ϕ , δ , ψ)| 2 β + n (θ + n , ϕ + )/β (θ , ϕ ), (12.133)
where |C + n | 2 for a reflected order of the number n in conical diffraction is expressed by

|C + n | 2 = ε + ε v |E + n | 2 + µ + µ v |B + n | 2 .
As mentioned in Sec. 12.3.2, the balance requirement (12.132) is one of the most important accuracy criteria based on a single computation generalized in the lossy case by the explicit computation of A from (12.73). The sum R + A is actually the energy balance for an absorbing grating in conical diffraction, including that having rough grooves, and the extent to which it approaches unity is a measure of the accuracy of a calculation.

For λ /d ≪ 1 the discrete order efficiencies is an approximation of the differential reflection coefficient (DRC) ς (analogous of a bistatic scattering coefficient [12.38]) for a continuum of scattered angles so that

∑ β + n >0 η + n = ∫ π/2 -π/2 ς (θ + n )dθ + n . ( 12
ς (θ + n ) = d cos 2 θ + n |C + n | 2 cos ϕ λ cos θ . ( 12.139)
The general case of 2D rough surfaces may be considered in a similar way. It can be done, for example, by expressing the solution of the 3D Maxwell equations for bigratings through solutions of the 2D Helmholtz equation for classical gratings working in conical diffraction, an approach which may be resorted to in some important cases (see Sec. 12.8.3).

Scattering intensity of rough gratings in a dispersive plane

For accounting random roughness rigorously in our codes, we use the model in which the randomly rough surface is represented by a grating of large period d. This period may contain a few or a large number of random asperities and/or a discrete number of periodic grooves. So the program deals with a structure that is a grating from a mathematical point of view but that can model a randomly rough surface of a grating or a mirror. If the groove spacing becomes small compared with the correlation length ξ of the random asperities, then the discrete dimension scaling can be applied to such a rough grating and the diffraction is investigated on the equivalent surface structure in proportionally higher diffraction orders. Moreover, if the width of the asperities has the same order of magnitude as the wavelength of incident light, the number of diffraction order is large, and the continuous speckle of the randomly rough surface is simulated by the discrete speckle of the grating, as has been demonstrated above.

In order to compute the scattering properties of a random rough surface using electromagnetic solvers (Penetrating or Separating) and a Monte Carlo procedure, an ensemble of surface realizations must be generated. There are several ways to generate a statistically stationary random surface [12.48]. The most common approach consists in generating surface profiles by the following technique. A sequence of random numbers (∼ 10 5 ) with Normal statistics, zero mean, and variance (rms roughness) σ = 1 is constructed from another random series directly generated by a computer. Then the former sequence is scaled in order to obtain a desired σ and, further, correlation with the Gaussian is performed in order to obtain a profile with a Gaussian correlation function. This is known as the spectral method ( [12.38]) and is used in PCGrate codes.

The boundaries of such randomized grating or mirror have both periodical and random roughness components and some averaging of random samples (from a few up to a several hundred) is required to obtain the exact scattered field intensities (see Example 10 in Sec. 12.9). In some conditions, fortunately not in x-ray-EUV, for instance when surface waves can propagate (like polaritons for metallic surfaces in the TM polarization), very big numbers of discritization points and propagating and evanescent orders (about a few thousand or even more) must be taken into account. Rigorous computation of the field scattered by random rough surfaces is a problem of daunting complexity in the area of electromagnetism and optics even for modern computers because of the very small wavelength-to-period and small wavelength-to-height ratios. It is especially true for x-ray-EUV grating and mirror applications. Therefore, the hardest diffraction problems may require large amounts of computer memory and, especially, high speed of computations.

Scattering intensity of rough gratings in a non-dispersive plane

The IMs, which have been developed in the frame of electromagnetic theory, permit application of optical grating methods to analysis of specular and diffuse x-ray-EUV scattering from rough gratings and mirrors using Monte Carlo calculus. The question of the closeness of results for 1D and 2D surfaces is of interest of this Section, since numerical methods for 1D surfaces are well established and efficient, and widely used for surfaces with 2D roughness [12.39]. The derivations of the boundary integral equations using potential operators as well as some details of their numerical implementation were described in previous Sections. An important case of bi-periodic gratings and 2D rough surfaces may be considered in a way by expressing the solution of the 3D Maxwell equations through solutions of the 2D general Helmholtz equations in conical diffraction, an approach which may be resorted to in short waves and shallow surface using the equivalence rule derived in App. D.

The effect of roughness on the 2D DRC can be exactly taken into account with model in which an uneven surface is represented by a bigrating with large periods of d x,z in perpendicular planes, which include appropriate numbers of random asperities with correlation lengths of ξ x,z . We analyze a complex structure which, while being the bigrating from a mathematical viewpoint, is actually the rough surface for d x,z ≫ ξ x,z . If ξ x,z ∼ λ and the number of orders is large, the continuous angular distribution of the energy reflected from randomly rough boundaries can be described by a discrete distribution η mn in orders (m, n) of a bigrating, similar to (12.134) for classical gratings. A study of the scattering intensity starts with obtaining statistical realizations of profile boundaries of the structure to be analyzed, after which one calculates the DRC for each realization, to end with the DRC averaged out over all realizations to obtain a mean DRC. By selecting large enough samples and numbers of sampling points, one comes eventually to properly averaged properties of the rough surface; however, this approach does not involve approximations, including averaging by the Monte Carlo method.

The equivalence rule between 2D and 1D grating efficiencies

A general approach to find efficiencies of bigratings and mean DRCs of rough 2D surfaces which permits one to use exact integral equations, rigorous (extended) boundary conditions, and radiation conditions leads to tedious calculus even in a case of perfectly conductive surfaces [12.49]. However, a great deal of simplification of the given boundary-problem can be achieved for shallow gratings and randomly-rough surfaces if we use the Rayleigh hypothesis together with the small-amplitude perturbation technique. Implementations of such a method, in which the reduced Rayleigh equations for reflection from such structure are solved in the form of expansions of the amplitudes of the p-and s-polarized components of the scattered field in powers of the surface profile function through terms, up to the third order, were proposed in several papers (see, e.g., Ref. 12.50 and references therein). In the present work, the authors use the perturbative analysis results only in order to derive an approximate connection rule between the efficiency of a shallow bigrating and efficiencies of two classical gratings with grooves rotating on 90 deg. The efficiency itself of a classical grating working in conical diffraction is defined rigorously using the boundary integral equation method, as it is prescribed in previous Sections.

The equivalence rule can be formulated as the following (see (12.163) of App. D)

η mn = η m η n r F , m ∨ n = 0, h x,z /d x,z < 1, (12.140)
where η m and η n are classical grating efficiencies obtained in conical diffraction, h x,z -profile heights in perpendicular planes, r F -the Fresnel factor of a 2D surface. It is worth noting that η m and η n in this equivalence rule should be computed with preservation of incidence and polarization angles of both gratings in the absolute coordinate system. Thus, using (12.140) the efficiency η mn of bigratings can be easily expressed in terms of the product of the efficiencies of two respective classical gratings considered in perpendicular dispersive planes and working in conical mounts at any polarization state. Equation (12.140) was derived in Ref. 12.37 for the normal incidence of linearly-polarized light on a simple boundary-profile bigrating. The equivalence rule described above is very similar to the impulse approximation result of the atomic scattering theory and can include multiple scattering in each perpendicular direction but always excludes cross-correlation components.

The derived connection equation is approximate and valid for shallow periodic surfaces of the type considered. However, this equivalence rule was checked successfully against various numerical examples, including non-shallow bigratings working at different wavelengthto-period ratios [12.37, 12.51]. It was found that it gives accurate results under the following assumptions: (a) h x,z d x,z and (b) λ d x,z . However, for non-deterministic surface profiles working in short waves, some modification of these conclusions is required. As follows from the known results obtained from analytic and asymptotic expressions valid for x-rays (see, e.g., Refs. 12.6, 12.40), (12.140) gives high-accuracy solutions for shallow rough 0D (i.e. rows of atoms with displacements), 1D, and 2D surfaces if the following conditions are fulfilled: (c) cos θ ′ h x,z ≪ d x,z and (d) λ ≪ d x,z cos ϕ , where θ ′ is an incidence angle on the surface. In case of x-ray-EUV ranges refractive indices of materials are close to the vacuum refractive index and h x,z can be large enough for grazing incidence. Thus (a) and (c) are close due to the nature of the perturbative development. However, (d) extends the range of the validation of (12.140) significantly, i.e. to the whole short-wave optical range because of the absence of optical resonances (i.e. due to plasmons, polaritons, waveguide resonances, etc) in x-rays and EUV.

Examples of numerical results

The described theoretical and numerical approaches for the calculation of far-zone fields and polarization properties of diffraction gratings are well suited to various types of optical grat-ing analysis. In this Section, we are going to analyze numerically examples of diverse grating diffraction problems. The results presented demonstrate the impact of diffraction and polarization incident angles, boundary shapes and layer refractive indices on diffraction and absorption in periodical structures. 12.2: Diffraction efficiencies (η ± ), diffraction θ ± , ϕ ± ) and polarization(δ ± , ψ ± ) angles of a dielectric lamellar grating a In this Section, we present several numerical experiments taken from well-known spectroscopic applications of gratings working in various mounts and polarization states at different wavelengths. More specifically, they are: the typical dielectric and metallic lamellar gratings illuminated in conical diffraction; the typical dielectric sine grating working in off-plane mounts; the typical metallic echelette gratings illuminated in conical diffraction; the anomalously absorbing Ag shallow-sine grating working in off-plane mounts in the visible; the photonic crystals with Au nanorods of various cross-sections illuminated at normal incidence in the visiblenear-infrared; the photonic crystal with dielectric circular nanorods working in different mounts in the near-and mid-infrared; the Al echelle grating protected by a thin layer of MgF 2 and illuminated in conical diffraction in the vacuum ultraviolet (VUV); the Au blaze grating working in grazing-incidence off-plane mounts in soft x-rays; the minimally-absorbing Mo/B 4 C multilayer blaze grating illuminated in grazing conical diffraction in soft x-rays; the flight Mo/Si multilayer trapezoidal grating working in the near-normal-incidence EUV and with random roughness accounting. The numerical examples of calculation results described in this Section were calculated using a few commercial and non-commercial IM-based codes.

DO b θ ± , • ϕ ± , • η ± ,% δ ± , • ψ ± , • R -2

Efficiencies and polarization angles of lamellar gratings

The efficiency results of reflection orders of the present IM for a typical conducting lamellar grating with the ridge width c and depth 2H working in a conical mount are demonstrated in Table 12.1. The grating and light parameters are as follows: c/d = 0.5, 2H/d = 1, ε + = 1, ε -= (-24.99, 1), µ ± = 1, λ /d = 0.5, θ = 22.208 • , ϕ = 20.705 • , δ = 45 • , and ψ = 0. We used 400 discretization points, mesh grading and the discretization of V + J -to calculate this example that allocates 188 MByte memory. The energy balance error calculated from (12.74) is ∼ 10 -6 . The average time taken up by the example on a workstation with two Quad-Core Intel Xeon 2.66 GHz processors, 8 MB L2 Cache, 1333 MHz FSB and 16 GB RAM is ∼ 1.5 s when operating on Linux Ubuntu 12.04 LTS 64-bit or Windows Vista Ultimate 64-bit and employing eightfold paralleling.

Table 12.3: Diffraction efficiencies (η ± ) and diffraction (θ ± , ϕ ± ) and polarization (δ ± , ψ ± ) angles of a dielectric sine grating for B z = 0 a

DO b θ ± , • ϕ ± , • η ± ,% δ ± , • ψ ± , • R -3 -
a 2H/d = 0.3, ε + = 1, ε -= 4, µ ± = 1, λ /d = 0.5, θ = 60 • , ϕ = 15 • , δ = 81.501 • , ψ = 0. b Diffraction order
In Table 12.2, the efficiency data of reflection and transmission orders for a similar dielectric lamellar grating in a conical mount are presented. The grating and light parameters are as follows: c/d = 0.5, 2H/d = 0.5, ε + = 1, ε -= 2.25, µ ± = 1, λ /d = 0.5, θ = 35.264 • , ϕ = 30 • , δ = 45 • , and ψ = 90 • . We used N = 400, mesh grading and the discretization of V + J -to calculate this example that allocates 188 MByte memory. The energy balance error calculated from (12.74) is ∼ 10 -5 . The average time taken up by the example is ∼ 1.5 s when operating on the aforementioned workstation and operating system. The efficiencies and polarization angles obtained in this and two next Subsections for transmission and reflection gratings working in conical diffraction can be compared with those obtained by the use of other rigorous methods and codes [12.7, 12.8].

Efficiencies and polarization angles of dielectric sine grating

In Tables 12.3 and 12.4, the efficiency results of the IM for a typical dielectric sine grating working in a conical mount are presented. The grating and light parameters are as follows:

Table 12.4: Diffraction efficiencies (η ± ) and diffraction (θ ± , ϕ ± ) and polarization (δ ± , ψ ± ) angles of a dielectric sine grating for E z = 0 a 12.3, the incident polarization angles are δ = 81.501 • and ψ = 0, for Table 12.4-δ = 8.499 • , ψ = 180 • . We used 100 discretization points and the numerical differentiation of V + to calculate these examples which allocate 10 MByte of RAM. The energy balance error calculated from (12.74) is about 10 -5 for both components of the polarization incident radiation. The average computation time taken up by an example on the aforementioned workstation and operating system is ∼ 0.1 s.

DO b θ ± , • ϕ ± , • η ± ,% δ ± , • ψ ± , • R -3 -
a 2H/d = 0.3, ε + = 1, ε -= 4, µ ± = 1, λ /d = 0.5, θ = 60 • , ϕ = 15 • , δ = 8.499 • , ψ = 180 • . b Diffraction order 2H/d = 0.3, ε + = 1, ε -= 4, µ ± = 1, λ /d = 0.5, θ = 60 • , ϕ = 15 • . For Table

Efficiencies and polarization angles of metallic echelette grating

The numerical results for a typical metallic echelette grating with blaze angle ζ and an apex angle of 90 • (see Fig. 12.2) working in a conical mount are demonstrated in Tables 12.5 and 12.6 for the two basic states of the incident polarization: δ = 0, ψ = 180 • or δ = 90 • , ψ = 0. The grating and light parameters are as follows: ζ = 30 • , ε + = 1, ε -= (-45, 28), µ ± = 1, λ /d = 0.5, θ = 0, ϕ = 40 • , and ψ = 0. One has used N = 800, mesh scaling near edges and the differentiation of V + to calculate these examples allocating 196 MByte of RAM. The average energy balance error calculated from (12.74) is ∼ 10 -5 for both polarization states of the incident radiation. The average computation time taken up by two values of the polarization angle on the aforementioned workstation and operating system is ∼ 3 s.

Anomalously absorbing Ag shallow-sine grating in the visible

Resonance and non-resonance anomalies differing in their nature can be effectively explored in high conductive gratings, such as: surface plasmon excitations, Bragg and Brewster conditions, groove shape features, etc. Because the s and p modes in conical diffraction are coupled through the boundary conditions, the associated problems are more general, and gratings act as perfect absorbers and local-field enhancers.

δ = 0, ψ = 180 •a DO b θ + , • ϕ + , • η + ,% δ + , • ψ + , • R -1 -
a ζ = 30 • , ε + = 1, ε -= (-45, 28), µ ± = 1, λ /d = 0.5, θ = 0, ϕ = 40 • . b Diffraction order
δ = 90 • , ψ = 0 a DO b θ + , • ϕ + , • η + ,% δ + , • ψ + , • R -1 -
a ζ = 30 • , ε + = 1, ε -= (-45, 28), µ ± = 1, λ /d = 0.5, θ = 0, ϕ = 40 • . b Diffraction order
In Fig. 12.17, the absorption of the Ag sinusoidal grating with d = 2.2 µm and 2H = 100 nm is calculated for the δ = 90 • , ψ = 0 or δ = 0, ψ = 180 • polarized incidence light with λ = 663 nm as a function of θ for ϕ = 0 (classical, TE and TM) or ϕ = 50 • (conical). The refractive indices of Ag were taken from Ref. 12.36 (µ ± = 1). For in-plane diffraction, anomalous absorption exists only for the TM polarization, while for conical diffraction both components are absorbed but in smaller amounts.

Note that we used the variant of discretization of H + V -to calculate these examples. The calculated problem allocates 10 MByte of RAM using N = 100. The energy balance error calculated from (12.74) is about 10 -6 for both components of the polarization of incident radiation. The average computation time taken up by the example on the aforementioned workstation and operating system is less than 0.1 s per calculation point.

Photonic crystals with Au nanorods in the visible-near-IR

In this Subsection, we are going to analyze numerically the optical response (reflection and absorbtion) of photonic crystal slabs supporting polariton-plasmon propagation with different cross sections of nanowires invariant with respect to the z axis. The essential physics of the formation of localized plasmon polariton modes in metallic nanowire arrays is described in Chapter 1. The vital role of the absorption, slab cross-section shape, and filling ratio of photonic crystals in the visible and near infrared regions is demonstrated in this Subsection. The model contains M -2 (see Fig. 12.3) identical gratings with closed boundaries (inclusions) of simple cross sections displaced vertically (by h m ) and horizontally (by f m ) relative to one another and embedded in a homogeneous medium with dielectric permittivity ε 1 and magnetic susceptibility µ 1 . We deal here only with materials with µ m = 1, m = 0, . . . M, although the model is applicable to other cases as well, including metamaterials [12.18]. The dependence of the dielectric permittivity ε m , m = 2, . . . M -1 of the material of nanorods on the incident photon frequency is assumed to be known. The lower medium and the upper medium are likewise assigned pairs of material constants, but one may conceive of more complicated cases of multilayer structures as well. The model also allows arbitrary incidence of, in the general case, elliptically polarized radiation on photonic crystals, which is prescribed by two angles of incidence and two angles of polarization. Figure 12.18 displays for comparison theoretical spectra of energy reflected from, and absorbed by, a photonic crystal with Au nanowires of circular, square, rectangular, and triangular cross sections of the same area and with M = 3 studied in the 1-3-eV photon energy range (visible and near infrared). In this and similar subsequent examples, we consider the TM-polarized (θ = ϕ = δ = 0, ψ = 180 • ) light normally falling on Au nanowires embedded in a SiO 2 matrix with d = 200 nm, ε 0 = ε 1 = ε 3 = 2.13, and refractive indices of Au taken from Ref. 12.36. The orientation of the rods having edges is chosen in such a way that light normally falls on one side of the rods. The a × b dimensions of the rectangular rods selected for this example are 50 × 25 nm 2 or 25 × 50 nm 2 and the width of the squares or triangles and diameter of the circles were chosen to obtain equal cross sectional area S = 1250 nm 2 . As seen from Fig. 12.18, reflection and, particularly, absorption spectra exhibit a strong difference near the plasmon-polariton anomaly among the five shapes of the nanowire cross section chosen. These differences amount to several hundred percent for the rectangles because of their different width-to-height ratio (two and a half) compared with the square or the circle (one) and the equilateral triangle (0.866). One observes also a noticeable difference in the positions of the absorption and reflection maxima among different grating profiles. Thus, the simple effective medium theory cannot be applied to design and analysis of such photonic crystals, even for a small filling ratio [12.13].

Figure 12.19 presents energy spectra similar to those displayed in Fig. 12.18 but for S four times that of the preceding example. In this case, a × b = 100 × 50 nm 2 or 50 × 100 nm 2 . We readily see that the differences in the reflection and absorption spectra among gratings of different profiles increase with increasing filling ratio and are observed now not only close to the plasmon resonances. Near the resonances, they amount to a few dozen percent of energy. The absorption spectra of the triangular-shaped nanowires have an interesting step-like function behaviour, which is not the case for absorption spectra of nanowires of other rod shapes.

Only 50 discretization points, mesh grading, Hankel kernel functions for inclusions and discretization of H + V -have been used to compute these examples which allocate ∼ 0.1 MByte memory. The relative error calculated from the energy balance for absorption gratings is ∼ 10 -4 . The average time taken up by one point on the aforementioned workstation and operating system is less than 0.1 s. 12.9.6 Lossless photonic crystal with circular rods in the near-and mid-IR

In this example, we consider numerically some diffraction properties of non-absorbing photonic crystals with dielectric rods. The influence of the geometry and number of crystal layers, the shape of rods, the filling ratio, the index of refraction of materials and the polarization and diffraction angles of light can be investigated for this type of photonic crystals. The role of the filling ratio, refractive index and polarization was demonstrated for the classical diffraction [12.12, 12.31]. Here we demonstrate, as an example of possibilities of developed software, the vital role of the filling ratio and polarization for conical diffraction. In Fig. 12.21 for the off-plane diffraction ϕ = 30 • and this is an additional parameter compared with the classical diffraction case.

For both in-plane and off-plane examples, there is a very different behavior in diffraction properties for TE and TM polarization components of the incident radiation, especially for big filling ratios. Compared with respective curves obtained in Figs. 12.20 and 12.21, it emerges that for s-polarized light the centers of the conical diffraction gaps have shifted significantly to smaller wavelengths and the widths and depths of the gaps have decreased considerably. In contrast to this behavior, for p-polarized light the centers of the conical diffraction gaps compared with the in-plane ones have shifted a little bit in opposite directions and the widths and depths of these gaps have increased noticeably. The vital importance of the azimuthal angle ϕ , as well as the incidence polarization has become evident even for a small filling ratio (p = 0.125); however they are more important for a high filling ratio (p = 0.5). Thus, using the conical diffraction for dielectric photonic crystals gives additional control parameters which significantly affect Bragg diffraction and existing photonic band gaps.

Only N = 50 without mesh grading and with Hankel kernel functions for inclusions are required to compute these examples using discretization of H + V -which allocates ∼ 0.2 MB memory. The relative error calculated from the energy balance for non-absorption gratings is ∼ 10 -4 . The average time taken up by one point on the aforementioned workstation and operating system is less than ∼ 0.1 s. 12.9.7 Al echelle grating coated by MgF 2 in the VUV Echelle gratings or simple echelles working in high spectral orders near Littrow diffraction conditions at high angles θ are one of the most popular grating types; however, they are rather difficult for fabrication and efficiency computations, especially those with dielectric coatings. A thin oxide film on the Al grating surface may lead to degradation of its diffraction properties at wavelengths below 130-140 nm. To protect and even improve the echelles' reflectance surfaces, a thin dielectric coating with a thicknesses of a few dozen nm can be applied in the VUV range. The usual material is MgF 2 , but sometimes other dielectrics are used. At a certain thickness of the coating film, waveguide phenomena come out to affect the grating performance; as a result, the diffraction efficiency can either decrease or increase as compared to the nonoxidized bare grating. A non-conformal layer which is obtained by two adjacent non-parallel boundaries (having different vertical distances between) provides a new freedom in design, but the analysis of gratings becomes more complex. Furthermore, echelles are frequently used in conical diffraction, making it possible to separate beams in a non-dispersive plane [12.3].

Our example deals with an aluminium echelle with 316 grooves/mm, working blaze angle ζ 1 = 63.4 • (r-2, i.e. tan ζ 1 = 2), and apex angle 90 • . The grating works at the -47th order, wavelength λ = 120 nm and ϕ = 6.5 • . A protecting MgF 2 layer is applied. Other coating and light parameters are as follows: ε 0 = 1, ε 1 = ε MgF 2 = (2.643876, 0), ε 2 = ε Al = (-1.2353087, 0.0913816), µ m = 1, m = 0, . . . 2, δ = 0, and ψ = 180 • . We consider four variants of its thickness and shape including zero thickness for the bare Al grating. For coated gratings, the coating's upper boundary is sawtooth, with right angle at the top vertex situated h 0 = 30 nm above the grating's top vertex. Thus, the variants differ from each other by the coating's working angle, which is case and ζ = 62.9 • or ζ = 63.9 • -for two non-conformal cases. Such non-conformal models of coatings do not pretend to be the best description of real structures formed by sputtering, but are simple and possible; and they account for a deviation of coating direction from the Al substrate surface which leads to a tapered shape on both slopes of the triangular profile. Fig. 12.22 presents angular dependencies of the grating efficiency. The efficiency results for such echelles obtained using different IM-based codes in in-plane mounts are presented in Ref. 12.25. Fig. 12.22 shows that the conformal coating leads to a noticeable increment of efficiency in comparison with a bare case over the whole range of angles, by ∼ 20%. The non-conformal coating with ζ = 62.9 • increases the efficiencies by ∼ 10% compared to the bare grating. The geometry in this case is such that the working facet receives a thinner layer of MgF 2 , which narrows approaching the vertex; the non-working facet gets a fatty coating. In contrast, the non-conformal coating with working angle ζ = 63.9 • does not increase the efficiency at its maximum and leads to practically the same efficiency graph as for the bare Al grating case in the whole central angular range. The opposite impact of these non-conformal coatings working in classical diffraction is demonstrated in Ref. 12.25. Thus, the efficiency is very sensitive to the boundary vertical shift, to the deviation of a MgF 2 layer from conformal shape, and also to the off-plane deviation.

Computations in this example were carried out with N = 800 for the bare grating and with N = 1600 for the gratings with conformal and non-conformal coatings. One also has used mesh scaling near edges and the differentiation of V + to calculate these examples, allocating 1024 MByte of RAM for N = 1600. The relative error calculated by (12.128) from the energy balance for absorption gratings is ∼ 10 -4 . In case of piecewise linear profiles, many pairs of kernel function arguments can be obtained from each other by translations; corresponding kernel function values are equal. Hence, there is an effective way to check for given arguments, whether or not we already encountered a congruent pair and calculated the kernel function for it (see Sec. 12.6.2). This approach significantly reduces computational time for echelles and even more-in case of conformal layers, where the kernel function values calculated on the upper side of the layer can be reused on the lower side. Calculation for each point on Fig. 12.22 required between a few s (bare case) and several dozen s (conformal and non-conformal cases) on the aforementioned workstation and operating system.

Au off-plane-grazing-incidence blaze grating in soft x-rays

The conical diffraction mount in which the direction of incident light is confined to a plane parallel to the direction of the grooves has the unique property of maintaining high and sustained diffraction efficiency, which is very important in the x-ray-EUV range. Such gratings are utilized as dispersive elements in laboratory and space spectral instruments, time-delayed compensators or splitters and spectral purity filters for EUV lithography. Grazing-incidence off-plane gratings have been suggested for the International X-ray Observatory (IXO) [12.53]. Compared with gratings in the classical in-plane mount, x-ray gratings in the off-plane mount have the potential for superior resolution and efficiency for the IXO mission. The results of efficiency calculations for such a gold blazed soft x-ray grating in a conical mount using the perfect triangular groove profile with d = 200 nm are shown in Fig. 12 In Fig. 12.23, the numerical results of the IM presented for a finite boundary conductivity are compared with those based on the IM with the perfect boundary conductivity multiplied by Fresnel reflectances calculated with respect to the blaze facet. The incident beam in the computations based on the perfect conductivity model and classical diffraction (using the In-variance theorem (see in [12.37] and Ch. 4) was assumed to be 100% TE-polarized (B z = 0). The refractive indices of Au were derived from the compilation at [12.55].

Rigorous computations carried out by the methods presented show that for the finite grating model all the order efficiencies are not sensitive to a polarization state. For both basic polarization state of the incident radiation order efficiencies presented in Fig. 12.23 differ not more than a few tenths of a %. Contrary, calculations based on the perfectly conducting boundary model are very sensitive to the polarization state and sharp Rayleigh anomalies for the TM-polarized incident radiation (not shown) occur. As can be seen in Fig. 12.23, the agreement between the data obtained by the finite conductivity model and the perfect conductivity model is good when the TE-polarization is used for the perfect conductivity model. The same conclusions were derived for a similar grating problem in Ref. 12.7 using the real (measured) average groove profile for the efficiency computation.

We have used 800 discretization points, the numerical differentiation of V + and no mesh scaling to calculate the finite-conducting blaze-groove-profile example that allocates a space of 144 MByte. The energy balance error calculated from (12.74) is ∼ 10 -4 in the investigated wavelength range. The average computation time taken up by one wavelength on the aforementioned workstation and operating system is ∼ 2 s. The time of a computation using the perfect conductivity model for N = 200 is about eighty times shorter at the same computation accuracy.

Efficiency of order

Wavelength, nm #-2, perfect cond., TE #-1, perfect cond., TE # 0, perfect cond., TE #+1, perfect cond., TE #-2, finie cond.

#-1, finie cond. # 0, finie cond. #+1, finie cond. 12.9.9 W/B 4 C multilayer off-plane-grazing-incidence blaze grating in soft-x-rays

Multilayer coated blazed gratings with high groove density are the best candidates for use in high resolution EUV and soft x-ray spectrometry such as resonance inelastic x-ray spectroscopy. Theoretical and experimental analysis show that such a grating can be potentially optimized for high dispersion and spectral resolution in a desired high number diffraction order without significant loss of diffraction efficiency. In order to realize this potential, the grating should have a perfect triangular groove profile and its absorption should be minimized. The grazing-incidence conical-diffraction mounting in which the direction of incident light is confined to a plane parallel to the direction of the grooves has the unique property of maintaining a maximal level of diffraction efficiency due to an additional angular parameter. In this Subsection, we analyze the optical absorption of a blazed multilayer grating working in grazing conical diffraction in the soft x-ray range.

In Fig. 12.24, the absorption of the 10000 /mm blazed Si grating coated with 60 bi-layers of W/B 4 C is calculated for the polarized (δ = 90 • , ψ = 0) incidence radiation with λ = 1.3 nm and θ = 6 • as a function of the azimuthal angle ϕ . The grating has a triangular groove profile with the blaze angle of 6 • and antiblaze angle of 64.53 • and a conformal multilayer coating (see Sec. 12.9.7) with the thicknesses of W and B 4 C layers measured in respect to the working facet normal, 0.6006 nm and 2.4024 nm , respectively. The refractive indices of Si, W, and B 4 C were taken from [12.55]. Figure 12.24 displays for comparison theoretical absorption spectra of a Si mirror coated with the same multilayer and working in the same mount. As one can see in Fig. 12.24, for the defined polar angle the grating and mirror absorptions are close in the azimuthal angle range investigated. Grating absorption minima less than 70% can be obtained for the azimuthal angle of ∼ 77.2 • . Thus, almost the all reflected energy can be directed into diffraction orders without additional losses for the multilayer soft-x-ray grating absorption.

Only N = 400 was used to compute this grating example accounting 121 boundaries which allocates ∼ 60 MB of RAM. The relative error calculated from the energy balance using Eq. (12.128) is ∼ 10 -4 . The average time taken up by one point on the aforementioned workstation and operating system is ∼ 4 min. L. Goray The depth of all the boundary profiles of the multilayer grating was 6.0 nm, with side slopes of 35 • and equal top and groove widths, as derived from the AFM and efficiency measurements. Because polarization effects are small near normal incidence, the efficiencies are presented for the case of TM-polarized radiation (ϕ = δ = 0, ψ = 180 • ). To determine the absolute values of order efficiencies, a model of the two-period-randomized-trapezium grating describing the realistic boundary shape and roughness was applied. For a rigorous accounting of the random roughness impact on the efficiency, the model with 41 randomly rough borders of the period of ∼ 476.19 nm having 400 random sampling points on two trapezoidal grooves with the same Gaussian surface roughness height statistics and Gaussian autocorrelation function was applied (for random border generation on non-flat surface shapes, see [12.1]) . The rough boundary parameters are as follows: the Si-Mo interface rms roughness σ Si-Mo = 0.2 nm and the Mo-Si rms roughness σ Mo-Si = 0.85 nm. The lateral correlation length ξ = 5 nm was chosen from the detailed microscopic analysis and the growth model of typical Mo/Si layers obtained by using magnetron sputtering [12.47]. An assumption about the absence of a vertical correlation between the border random roughness components was applied in this model. Seven sets of 41 rough border profiles were generated to compute exact efficiencies of the FL1 multilayer grating. The Si protective capping layer of 2 nm was modeled by using 1.5-nm-thick To determine the absolute values of scattering light intensities between orders [12.5], a model of ten-period-randomized-trapezium grating allowing a fine-diffraction-angle discretization and describing the realistic boundary shape and roughness was applied. For a rigorous accounting of random roughnesses, the model with 41 randomly rough low-frequency borders of the period of ∼ 2381 nm having 800 random sampling points on 10 trapezoidal grooves with the described above rough boundary parameters was used. Some 105 sets of 41 non-correlated vertically border profiles were generated to compute exact scattering light intensities between orders of the FL1 multilayer grating. The same layer parameters as for the efficiency model (see above) were used. Refractive indices derived from the NIST data for Mo [12.57], from the CXRO data-for Si [12.53], and from the Palik data-for SiO 2 were used for efficiency calculations in the whole wavelength range.

Convergence and accuracy of the efficiency results of the randomly-rough 41-boundary grating were investigated using the Penetrating solver, Gauss computation algorithm and finite type of low border conductivity. All the accelerating convergence options were switched on in PCGrate-SX v. 6.5 (see Sec. 12.4.3). The linear type of refractive index data interpolation was chosen. A high rate of convergence of the results was observed for the developed grating efficiency model [12.1]. Only several sets of 41 rough border profiles and the medium number of discretization points per boundary are enough to compute exact efficiencies in all orders of interest. The differences between principal order efficiencies obtained with seven boundary sets with low (N = 600), medium (N = 800), and high (N = 1000) accuracy are about a few percents for all orders under study. The differences between efficiencies obtained with three, five, and seven statistical boundary sets (N = 800) are also about a few percent for all diffraction orders under study. For the final efficiency modeling (Fig. 12.25), N = 800 and seven random boundary sets are used. The total error for all points and ranges derived from the energy balance (12.128)

Appendix D: Derivation of the general connection rule between 2D and 1D gratings

We seek for a perturbative development of the reflection operator R in powers of the heights h (i)

x and h ( j) z of a bi-periodic surface either conductive or dielectric that is the sum of the two Fourier series:

h(x, z) = h x + h z = ∑ i h (i)
x sin(2πxi/d x + τ Such a representation of h(x, z) is typical for real 2D periodic or random surfaces obtained, e.g., as a linear response of a photoresist to light with two separate exposures in perpendicular planes or by polishing using a linear tool. Note that the 2D Fourier transformation of h(x, z) is also the sum of two 1D Fourier transforms of h x and h z . We suppose also that the bigrating works under arbitrary incidence and polarization states of a plane monochromatic wave and the respective single-periodic gratings work in conical diffraction. Suppose for simplicity h x and h z are even functions, which is true for many ergodic stationary processes. So, replacing h x or h z by -h x or -h z does not change the diffraction pattern in the far-field zone. We will study the perturbative expansion of the reflected efficiency η as a function of the surface heights (h (i)

x ) 2 and (h ( j) z ) 2 . Using the perturbative expansion of R, the terms of η which contain an expression such as (h η mn = η m η n r F + o(h 4 ), m ∨ n = 0, h x,z /d x,z < 1, (12.163) where η m and η n are 1D grating efficiencies obtained in conical diffraction, r F -the Fresnel factor of a 2D surface. It is worth noting that η m and η n in this equivalence rule should be computed with preservation of incidence and polarization angles of both gratings in the absolute coordinate system.

Figure 12 . 1 :

 121 Figure 12.1: Schematic conical diffraction by a grating.

Figure 12 . 2 :

 122 Figure 12.2: Cross section of a simple grating of period d with incidence direction k k k, incidence angle θ and conical angle ϕ .

Figure 12 . 3 :

 123 Figure 12.3: Cross section of a multilayer grating with inclusions and separating boundaries.

Figure 12 . 4 :

 124 Figure 12.4: Cross section of a multilayer grating with penetrating boundaries.

Figure 12 . 5 :

 125 Figure 12.5: The computation time for the lamellar grating example described in Table12.1.

Figure 12 . 6 :

 126 Figure 12.6: The computation time for a typical x-ray grating efficiency problem with repeated pairs of layers vs. number of identical boundaries.

Figure 12 . 7 :Figure 12 . 8 :Figure 12 . 9 :Figure 12 . 10 :Figure 12 . 11 :

 12712812912101211 Figure 12.7: Energy balance error with the standard IM using acceleration convergence terms for the problem of diffraction on a plane transmission interface (normal incidence in vacuum with the lower medium refractive index ν 1 = 1.5) vs. N for λ /d = 1.

LFigure 12 Figure 12 . 13 :

 121213 Figure 12.12: Reflected -1-order efficiency of an Au sinusoidal 300 grooves/mm grating with a depth of 25 nm illuminated by non-polarized radiation with λ = 4.4 nm incident at θ = 87.35 • , plotted vs. N for the standard IM or the MIM.

Figure 12 . 14 :

 1214 Figure 12.14: Energy balance and transmitted energy with the MIM, plotted vs. g at N = 100 for the same diffraction problem as in Fig. 12.7 but for λ /d = 0.01.

LFigure 12 . 15 :

 1215 Figure 12.15: Energy balance of an Au sinusoidal 300 grooves/mm grating with a depth of 10.5 nm illuminated by non-polarized radiation with λ = 0.834 nm incident at θ = 88.65 • , plotted vs. g for N = 100 and the standard IM or the MIM.

Figure 12 . 16 :

 1216 Figure 12.16: Energy balance vs. g used for the same diffraction problem as in Fig. 12.15, but for N = 400.

  a c/d = 0.5, 2H/d = 1, ε + = 1, ε -= (-24.99, 1), µ ± = 1, λ /d = 0.5, θ = 22.208 • , ϕ = 20.705 • , δ = 45 • , ψ = 0. b Diffraction orderTable

  a c/d = 0.5, 2H/d = 0.5, ε + = 1, ε -= 2.25, µ ± = 1, λ /d = 0.5, θ = 35.264 • , ϕ = 30 • , δ = 45 • , ψ = 90 • . b Diffraction order

Ligure 12 . 17 :

 1217 . Goray and G. Schmidt: Boundary Integral Equation Methods for Conical Diffraction ... Absorption of an Ag sinusoidal grating with d = 2.2µm and a depth of 100 nm working in classical (ϕ = 0) or conical (ϕ = 50 • ) diffraction, plotted vs. θ for λ = 663 nm and δ = 90 • , ψ = 0 or δ = 0, ψ = 180 • .

Figure 12 . 18 :

 1218 Photon energy, eV

Figures 12 .

 12 Photon energy, eV

Figure 12 . 19 :

 1219 Figure 12.19: The same as in Figure12.18, but for S = 5000 nm 2 .

Figure 12 . 20 :

 1220 Figure 12.20: Calculated transmission spectra of 1 µm-period gratings with 15 layers of dielectric circular rods with ε = 4 and different filling ratios p embedded in vacuum with hexagonal structure, plotted vs. λ for θ = 0, ϕ = 0 and different polarization angles (classical diffraction).

Figure 12 . 21 :

 1221 Figure 12.21: The same as in Figure 12.20, but for ϕ = 30 • (conical diffraction).

  igure 12.22: Efficiency in -47 order of A1 316 grooves/mm echelles with blaze angle ζ 1 = 63.4 • working in conical diffraction at λ = 120 nm and ϕ = 8 • : bare (ζ = 0), or with MgF 2 coating having upper sawtooth boundary with vertical displacement h 0 = 30 nm and working angle ζ = 63.4 • (conformal case), or ζ = 62.9 • (non-conformal case), or ζ = 63.9 • (non-conformal case), plotted vs. θ .

. 23 .

 23 The design blaze angle ζ is 7.5 • and the technique anti-blaze angle is 64.53 • [12.54]. Remaining grating and light parameters are as follows: µ ± = 1, θ = 0, ϕ = 88 • , and δ = 90 • and ψ = 0 or δ = 0 and ψ = 180 • .

Figure 12 . 23 :

 1223 Figure 12.23: Diffraction efficiencies of an Au triangular-groove-profile grating with d = 200 nm, ζ = 7.5 • , µ ± = 1 and for the incident wave with θ = 0, ϕ = 88 • and δ = 90 • , ψ = 0 or δ = 0, ψ = 180 • , plotted vs. λ .

Figure 12 . 24 :

 1224 Figure 12.24: Absorption of structures with 60 W/B 4 C bilayers on Si for the polarized (δ = 90 • , ψ = 0) grazing incidence x-ray radiation with λ = 1.3 nm and θ = 6 • vs. ϕ . 1-mirror; 2--blazed grating with 10000 grooves/mm and ζ = 6 • .

Figure 12 . 25 :

 1225 Figure 12.25: Calculated TM efficiencies of orders of a 4200 grooves/mm rough trapezoidal grating with 20 Mo/Si bilayers on Si operating at θ = 6.5 • vs. λ .

Figure 12 . 26 :

 1226 Figure 12.26: Calculated TM scattering intensities of a 4200 grooves/mm rough trapezoidal grating with 20 Mo/Si bilayers on Si operating at θ = 6.5 • vs. λ .

  will be denoted R kl :η = R 00 + R 01 + R 10 + R 11 + R 02 + R 20 + ...Using the quasi-periodicity property of R and Taylor expansion of scattered field amplitudes in powers of the surface profile heights (e.g, see Eq. 53 of Ref. 12.50), η mn , η m , and η n can be expressed in the following forms:

η

  mn -o(h 6 ) = δ mn a 00 + ∑ δ m,n is the Kronecker delta. From (12.159)-(12.161) we choose from the physical point of view one of the two possible expressions for η mn through η m and η n :η mno(h 6 ) =η m η n a 00(12.162) one can formulate the equivalence rule:

  .63) and similar expressions for the boundary integrals involving B z .Note that ε + and µ + are nonzero real numbers, and let ε -and µ -also be real. Taking the imaginary part of Eqs. (12.61) and (12.62) one gets

  -y 2 /4t ∑

								n∈Z	e -α 2 n t e inx dt √ πt	.
	Poisson's summation formula gives
	∑ n∈Z	e -(α+n) 2 t e inx =	√	π t	m∈Z e -iαx-x 2 /4t ∑	e -π 2 m 2 /t e πmx/t e 2πimα ,
	which leads to
	Ψ e (x, y) =	e -iαx 4π ∑ m∈Z	e 2πimα	∫ a 2 0	e k 2 t e -((x-2πm) 2 +y 2 )/4t dt t	.	(12.90)

  To avoid numerical overflow problems, which may occur if |y|/a is large, we use the relation

	w (	aβ n -i	|y| 2a	)	= 2 e y 2 /4a 2	e -a 2 (k 2 -α 2 n ) e i|y|β n -w (	-aβ n + i	|y| 2a	)	(12.95)
	obtained from (12.93), which gives
	i e -y 2 /4a 2 e a 2 k 2 8π	e -a 2 α 2 n β n	(	w (	aβ n -i	|y| 2a	)	+ w (	-aβ n + i	|y| 2a

12.94)

From

(12.93)

it can be seen that | w(z)| = O(e (Im z) 2 -(Re z) 2 ) if Im z < -| Re z|.

  + xn = k sin θ + n cos ϕ and k x = k sin θ cos ϕ we know the derivative

	From the grating equation for conical diffraction in the form
			k + xn = k x +	2πn d	,	(12.135)
	where k dk + xn dn	=	2π d	.	(12.136)
	Then, for large enough N, |n| ≤ N and accounting dn = 1 one can derive
	∑ β + n >0	dn = ∑ β + n >0	=	d 2π	∫ π/2 -π/2	k cos θ + n cos ϕ dθ + n .	(12.137)
	From (12.133), (12.132), and (12.137) we have
	∑ β + n >0	|C + n | 2 β + n /β =	d 2π	∫ π/2 -π/2	k cos 2 θ + n |C + n | 2 cos ϕ cos θ	dθ + n .	(12.138)
	Compare (12.134) and (12.138) we obtain the DRC for conical diffraction
								.134)

Table 12 .

 12 1: Diffraction efficiencies (η + ), diffraction (θ + , ϕ + ) and polarization (δ + , ψ + ) angles of a metallic lamellar grating a

	DO b θ + , •	ϕ + , •	η + ,%	δ + , • ψ + , •
	R -2 -43.715	-20.705 7.52	61.85 48.30
	R -1 -9.007	-20.705 13.25	15.79 -12.23
	R 0	22.208	-20.705 44.27	41.33 170.15
	R 1	65.852	-20.705 31.05	75.64 166.30

Table 12 . 5 :

 125 Diffraction efficiencies (η + ) and diffraction (θ + , ϕ + ) and polarization (δ + , ψ + ) angles of a metallic echelette grating for

Table 12 .

 12 6: Diffraction efficiencies (η + ) and diffraction (θ + , ϕ + ) and polarization (δ + , ψ + ) angles of a metallic echelette grating for

  and G. Schmidt: Boundary Integral Equation Methods for Conical Diffraction ... 12.7312.9.10 Flight Mo/Si multilayer rough lamellar grating in the EUVHere we present examples of the Mo/Si lamellar grating efficiency standardized for the Extreme-Ultraviolet Imaging Spectrometer (EIS) on the Hinode (former Solar-B) mission[12.56], the first implementation of a multilayer grating on a satellite instrument. We describe the performance of the flight FL1 4200 grooves/mm multilayer grating operating at θ = 6.5 • of the in-plane configuration in the wavelength region 17-21 nm. The efficiency was calculated by PCGrate-SX v.6.5 software using data of AFM measurements and was compared to the synchrotron efficiency measurements[12.1].

							#+1, meas.
		0.12					# 0, meas.
							#-1, meas.
							#+1, model
							# 0, model
							#-1, model
	Efficiency of order	0.04 0.08
		0
		17	18	19	20	21	22
				Wavelength, nm

 12.75was on the order of 10 -3 . The time taken up by one rigorous computation (one scanning point) for N = 800 on the aforementioned workstation and operating system is ∼ 45 min.

Convergence and accuracy of the scattering intensity results of a randomly-rough 41boundary grating were investigated using similar accuracy parameters to the efficiency computation. A medium rate of convergence of the light intensity results was observed for a wavelength of 19.25 nm and the above computation model. More than 100 sets of 41 rough border profiles and medium number of discretization points are enough to compute exact values of scattering light intensities between orders. The difference between scattered light intensities data obtained with seven boundary sets using medium (N = 1000) and high (N = 1200) accuracy is about 10 -5 for almost all diffraction angles. The differences between scattered light intensities obtained with different numbers of statistical boundary sets [START_REF] Goray | Numerical analysis of the efficiency of multilayer-coated gratings using integral method[END_REF]70,98,105) and N = 1000 for all diffraction (scattering) angles are shown in Ref. [START_REF] Maystre | Photonic crystal diffraction gratings[END_REF].1. For the final scattering intensity modeling (Fig. 12.26), N = 1000 and 105 random boundary sets were chosen. The total error for all points and ranges derived from the energy balance was on the order of 10 -5 . The time taken up by one computation for N = 1000 on the aforementioned workstation and operating system is about two hours.

Appendix A: Derivation of the recursive algorithm for Separating solver

In any of the strips {u j < y < d j-1 } the functions ε and µ take constant values and we introduce its wave number κ j by κ 2 j = ε µε 0 µ 0 sin 2 ϕ . As quasi-periodic solutions of the Helmholtz equation

in the strips {u j < y < d j-1 } between Σ j and Σ j-1 , j = 1, . . . , M -1, the functions E z , B z are smooth and bounded. Hence, for y ∈ (u j , d j-1 )

Assign to each profile Σ j a characteristic y-coordinate y j , for example y j = Y j (0) for a given parametrization (X j (t),Y j (t)) of the profile Σ j . Recall that y 0 > y 1 > . . . > y M-1 . Using the notation

the field in {u j < y < d j-1 } above Σ j is given by

) e iα n x , (12.142) (y-y j ) represent the diffracted waves. Thus, the coefficients in equations (12.142) and (12.143) are linked by the reflection and transmission matrices of the grating having only the interface Σ j .

For a compact notation, we introduce the infinite coefficient vectors

Then equations (12.141) can be written in the form

with the infinite diagonal matrix

Denoting by r j , t j the (infinite) reflection and transmission matrices of the grating with profile Σ j for illumination from above and by r ′ j , t ′ j the corresponding matrices for illumination of Σ j from below. This means that the incoming field with coefficient vector A j is diffracted by the simple grating with profile Σ j into the reflected field with coefficient vector r j A j and the transmitted field with coefficient vector t j A j . Analogously, illumination from below by a field with coefficient vector B j results in a reflected field characterized by r ′ j B j and a transmitted field with coefficient vector t ′ j B j . Hence, for any j = 1, . . . , M -2 the coefficient vectors are linked by the relations

Writing (12.11) in the form

we obtain (12.145) with j = 0, whereas for y < -H we derive from

Here we provide the formulas for solving the multi-profile problem to determine the vectors B B B 0 and A M-1 from given input A A A 0 and vanishing B M-1 . The idea is to look for a recursion for the operators R j , T j such that

By (12.146) we know that R M-1 = r M-1 , T M-1 = t M-1 . Furthermore, we have from (12.144) and (12.145)

L. Goray and G. Schmidt: Boundary Integral Equation Methods for Conical Diffraction ... 12.77 which gives

The last equation implies

which transforms (12.147) into

and hence

Finally, from

we derive

This leads to the following marching procedure:

12.11 Appendix B: Derivation of the recursive algorithm for Penetrating solver

The scheme is based on the ansatz

with certain 2 × 2 linear operator matrices A j and B j . Note first that the initial values (12.124) follow from (12.120) and the jump relation (12.26) for ∂ n S Γ M -1,M . Using (12.151) the transmission conditions (12.116) on Γ j for j = 1, . . . , M -1 can be written in the form The representation (12.119) and the jump relation (12.26) of the double layer potential D Γ j , j imply that

Hence (12.153) leads, in matrix notation, to the equation

) -

) , (12.154) which is equivalent to (12.116). Using the singular integral H ( j) (12.31)) we obtain the relation ((

) ( φ j-1

) ,

which is satisfied by

) ,

provided that Q j-1 is a solution of the operator equation (12.123). The equations (12.125) and (12.126) for A j-1 and B j-1 are derived from relations on the upper boundary Γ j-1 of G j . The representation (12.119) and condition (12.153) give

) , which by (12.151), (12.121) and using H

j-1 j ∂ t leads to (12.125). Now (12.126) follows from (12.23) and (12.119), since

imply that on Γ j-1

) .

Equation (12.127) follows from the relations

on the upper profile Γ 0 , which hold because E i z , B i z satisfy the Helmholtz equation (∆+ω 2 κ 2 0 )u = 0 and the outgoing wave condition in G - 0 = R 2 \ G 0 . Hence, the transmission conditions (12.115) are fulfilled if and only if

) -

) ,

i.e., if φ 0 , ψ 0 satisfy (12.127).

Remark 12.11.1 If the material in the bottom layer G M is a perfect conductor, then the zcomponents of E and B have to satisfy the boundary condition

In this case, it is easy to see that the relations (12.125) and (12.126)

and the initial values

) ,

) .

Hence, the densities {φ j , ψ j }, j = 0, . . . , M -2, are derived by the same scheme (12.121 -12.127).

12.12 Appendix C: Derivation of the absorption energy for multilayer gratings

As in Section 12.3.3 the application of Helmholtz equations and Green's formula in Ω H ∩ G 0 implies the relation

where (E z , B z ) is the solution of the conical diffraction problem, and

are the normal derivatives on Γ 0 of the z-components of the total fields in G 0 , i.e. the sum of the reflected and the incident fields. Setting the energy of the incident wave

the sum of reflection order efficiencies R (cf. (12.67)) fulfils

12.80

Gratings: Theory and Numeric Applications, Second Edition, 2014

Hence, for non-transparent grating we derive the energy conservation

If, otherwise, the material parameters ε M and µ M are real, then some part of the inci- dent field will be transmitted. Then similar considerations in the domain Ω H ∩ G M lead to the relations

where T is the sum of transmission order efficiencies of the multilayer grating (cf. (12.68)), and

z are the normal derivatives on Γ M-1 of the z-components of the transmitted fields in G M . In this case we derive the energy conservation

where the absorption A is given by the formula

) .

(12.157)

Using the jump conditions the obtained formulas for A can be easily transformed. For example, from (12.115) we know that on Γ 0

) .

Further, on Γ M-1

) .