A Plane wave polarization B S-matrices of corrugated gratings C Derivation of formulas describing the light diffraction on corrugated

Résumé

Il y a actuellement un vif intérêt pour des méthodes rigoureuses qui effectuent l'analyse électromagnétique des milieux diélectriques avec une distribution de permittivité diélectrique complexe. L'intérêt est motivé par des applications actuelles et futures dans la conception et la fabrication d'éléments optiques et optoélectroniques. Le niveau que les technologies de microstructuration ont maintenant atteint requiert des appels pour méthodes numériques rapides, économes en mémoire et rigoureuses capables de résoudre et d'optimiser des grandes parties de structures dont les caractéristiques représentent la fonction optique de la structure complète.

Bien que la majorité des problèmes de modélisation en microoptique sont non périodiques (par exemple, une section d'une couche diffusante d'OLED, la cellule d'un réticule microélectronique, une microlentille diffractive de haute NA), ils peuvent être efficacement résolus par la périodisation de la distribution de l'indice. Une nouvelle méthode numérique puissante pour la modélisation exacte de structures périodiques 2D est décrite avec toutes les fonctionnalités et les expressions nécessaires à son exécution. La puissance de cette méthode est dans sa forme spécifique unique qui permet d'appliquer rapidement des algorithmes numériques et, par conséquent, de diminuer de façon spectaculaire la complexité de calcul en comparaison avec les méthodes établies. La comparaison avec des solutions de référence a montré que, d'abord, la nouvelle méthode donne les mêmes résultats que celles-ci sur les structures de référence et, d'autre part, que le temps de calcul nécessaire et le recours en mémoire représentent une percée vers la résolution de grandes structures périodiques ou périodisées.

La méthode développée a été appliquée à analyser le problème de diffusion non périodique d'une couche diélectrique plan avec micro/nanoparticules sphériques. Une référence numérique proposée a démontré la possibilité d'obtenir environ 1% de précision. En outre, il a été développé un modèle numérique basé sur des matrices S pour la simulation des structures plane électroluminescentes. La validité de la méthode a été démontrée par comparaison avec les résultats expérimentaux. Enfin, les deux méthodes de calcul de la diffusion de la lumiére et de simulation des structures multicouches ont été groupeées, et une couche diffusante a été démontrée augmentait l'efficacité externe d'une OLED de quelques pour cents. i

Resume

There is presently a strong interest for rigorous methods that perform the electromagnetic analysis of dielectric media with complex dielectric permittivity distribution. The interest is motivated by both present and future applications in the design and manufacturing of optical elements and optoelectronic devices. The level that the microstructuring technologies have now reached calls for fast, memory sparing, and rigorous numerical methods capable of solving and optimizing large structure parts whose characteristics do represent the optical function of the whole structure.

Although the majority of modeling problems in microoptics are non-periodic (e.g., a section of an OLED extraction layer, the cell of a microelectronic reticle, a high NA diffractive microlens) they can be efficiently solved by periodizing the index distribution. A new powerful numerical method for the exact modeling of 2D periodic structures is described with all features and expressions needed to implement it. The power of this method is in its unique specific form which permits to apply fast numerical algorithms and, consequently, to decrease dramatically the calculation complexity in comparison with established methods. The comparison with reference solutions has shown that, first, the new method gives the same results as the latter on benchmark structures and, secondly, that the needed calculation time and memory resort represent a breakthrough towards solving larger periodic or periodized structures.

The developed method was applied to analyze nonperiodic scattering problem of a plane dielectric layer with spherical micro/nanoparticles. Proposed numerical benchmark demonstrated the possibility to get about 1% accuracy. In addition there was developed a numerical S-matrix based method for planar electroluminescent structures simulation. Validity of the method was demonstrated by comparison with experimental results. Finally both methods for the light scattering calculation and multilayer structures simulation were joined, and a scattering layer was demonstrated to increase an OLED external efficiency by several percent.

Résumé substantiel

Il y a actuellement un vif intérêt pour des méthodes rigoureuses qui effectuent l'analyse électromagnétique des milieux diélectriques avec une distribution de permittivité diélectrique complexe. L'intérêt est motivé par des applications actuelles et futures dans la conception et la fabrication d'éléments optiques et optoélectroniques. Le niveau que les technologies de microstructuration ont maintenant atteint requiert des appels pour méthodes numériques rapides, économes en mémoire et rigoureuses capables de résoudre et d'optimiser des grandes parties de structures dont les caractéristiques représentent la fonction optique de la structure complète.

Différents modèles décrivant la diffusion de la lumière sur les particules individuelles et de groupes de particules considèrent généralement les diffuseurs placés dans un milieu homogène et isotrope ou périodiquement reproduits à l'infini. Ces modèles et méthodes connexes font face à certaines difficultés cependant quand le volume de diffusion est infinie en deux dimensions et délimité par des interfaces planes dans la troisième dimension. Les exemples sont l'impact d'une couche de diffusion sur l'efficacité des OLEDs et des éléments diffractifs complexes de grande ouverture.

Les interfaces planes dans la zone de champ proche de particules diffusantes peuvent être prises en compte par différentes méthodes rigoureuses. La plupart d'entre elles, par exemple, les différences finies, éléments finis ou méthodes des équations intégrales de volume, sont d'une grande complexité numérique. Le meilleur choix pourrait être la méthode de la matrice T qui a été appliquée pour des particules près de la surface des géométries. Cependant, cette méthode nécessite des efforts supplémentaires de calcul de la matrice T pour chaque particule dans un ensemble de diffusion.

Un autre moyen de calcul de la dispersion dans des structures planes a été établi au moyen de méthodes de calcul de diffraction de lumière sur des réseaux. Principalement, ce sont des méthodes de Fourier, et, en particulier, la méthode Fourier modale (FMM). Récemment, la FMM a été appliqué au calcul de la diffusion des ondes électromagnétiques sur des objets 2D. Un avantage important de cette approche est que la forme géométrique de l'objet n'affecte pas la complexité de la méthode ni le temps de calcul. Cependant, la FMM a une complexité assez élevée égale à O(N 3 O ) avec N O étant le nombre d'harmoniques dans la transformée de Fourier de l'espace.

La présente thèse propose une méthodologie qui permet de résoudre exactement des grands systèmes, passe autre la limite O(N 3 O ), et diminue également la taille de mémoire requise. Pour ce faire, elle calcule un grand système 2D-périodique en un temps proportionnel à N O . Ceci est réalisé par le calcul d'une équation intégrale qui est réduit à un système d'équations linéaires dans la forme d'un produit de matrices bloc-diagonales et bloc-Toeplitz. Le système est résolue par des algorithmes de calcul connus comme la FFT et la GMRES.

Cette nouvelle méthode numérique puissante pour la modélisation exacte de structures périodiques 2D est décrite avec toutes les fonctionnalités et les expressions nécessaires à son exécution. La puissance de cette méthode est dans sa forme spécifique unique qui permet d'appliquer rapidement des algorithmes numériques et, par conséquent, à diminuer de façon spectaculaire la complexité de calcul en comparaison avec les méthodes établies. La comparaison avec des solutions de référence (données par la FFM et les méthodes de Rayleigh) a montré que, d'abord, la nouvelle méthode donne les mêmes résultats que celles-ci sur les structures de référence et, d'autre part, que le temps de calcul nécessaire et le recours de mémoire représentent une percée vers la résolution de plus grandes structures périodiques ou périodiséer.

La méthode développée a été appliquée à l'analyse d'une problème de diffusion non périodique d'une couche diélectrique plane avec micro/nanoparticules sphériques. Une référence numérique avec la solution de Mie a démontré la possibilité d'obtenir environ 1% de précision. En outre, il a été développé un modèle numérique basé sur les matrices S pour la simulation de structures planes électroluminescentes. La validité de la méthode a été démontrée par comparaison avec des résultats expérimentaux.

Avec l'application de la méthode à l'analyse des OLEDs avec couche diffusante, le modèle de propagation des ondes planes dans les OLEDs a été révisé et sa capacité à simuler rigoureusement toutes les propriétés électromagnétiques des structures a été démontrée. Les relations importantes exactes pour le flux de puissance et les pertes dans les couches OLED ont été données. L'applicabilité de ce modèle a été confirmée par une comparaison avec les données expérimentales obtenues par la mesure des propriétés optiques d'une OLED fabriquée verte.

Enfin, les deux méthodes de calcul de la diffusion de la lumière et de simulation des OLEDs ont été groupées, et une couche de diffusion a été démontrée augmentant l'efficacité externe d'une OLED de quelques pour cents.

Substantial resume

There is presently a strong interest for rigorous methods that perform the electromagnetic analysis of dielectric media with complex dielectric permittivity distribution. The interest is motivated by both present and future applications in the design and manufacturing of optical elements and optoelectronic devices. The level that the microstructuring technologies have now reached calls for fast, memory sparing, and rigorous numerical methods capable of solving and optimizing large structure parts whose characteristics do represent the optical function of the whole structure.

Various models describing the light scattering on single particles and groups of particles usually consider scatterers placed in a homogeneous isotropic medium or periodically continued to the infinity. These models and related methods face certain difficulties however when a scattering volume is infinite in two dimensions and bounded by plane interfaces in the third dimension. Examples are calculation of a scattering layer impact on the efficiency of photovoltaic devices and complex high-aperture diffractive elements.

Plane interfaces in the near field zone of scattering particles can be taken into account by different rigorous methods. Most of them, e.g., finite-difference, finite-element or volume integral equation methods are of a high numerical complexity analyzing multiparticle large aperture 3D scattering structures. The better choice could be the T-matrix method which was applied for particle-near-surface geometries. However, this method requires additional T-matrix calculation efforts for each particle in a scattering ensemble.

Another way to the scattering calculation in planar structures was established by means of methods for the light diffraction calculation on gratings. Primarily, these are Fourier methods, and, in particular, the Fourier modal method (FMM). Recently, the FMM has been applied for the calculation of the electromagnetic wave scattering on 2D objects. A prominent advantage of this approach is that scattering object geometry does not affect the method complexity and calculation time. However, the FMM itself exhibits quite high complexity equal to O(N 3 O ) with N O being the number of harmonics in the Fourier-space.

The thesis proposes a methodology which exactly solves large systems and breaks through the O(N 3 O ) limit, and also decreases the required memory size. It does so and calculates a large 2D-periodic system in a time proportional to N O . This is achieved via the analytical derivation of an adequately formulated integral equation which is reduced to v a system of linear equations in the form of a product of block-diagonal and block-Toeplitz matrices. The system is processed by known powerful calculation algorithms based on the fast Fourier transform (FFT) and the generalized minimal residual method (GMRES).

A new powerful numerical method for the exact modeling of 2D periodic structures is described with all features and expressions needed to implement it. The power of this method is in its unique specific form which permits to apply fast numerical algorithms and, consequently, to decrease dramatically the calculation complexity in comparison with established methods. The comparison with reference solutions (given by the FFM and Rayleigh methods) has shown that, first, the new method gives the same results as the latter on benchmark structures and, secondly, that the needed calculation time and memory resort represent a breakthrough towards solving larger periodic or periodized structures.

The developed method was applied to analyze nonperiodic scattering problem of a plane dielectric layer with spherical micro/nanoparticles. Proposed numerical benchmark with the Mie solution demonstrated the possibility to get about 1% accuracy. In addition there was developed a numerical S-matrix based method for planar electroluminescent structures simulation. Validity of the method was demonstrated by comparison with experimental results.

With a view of applying the method to the analysis of organic light-emitting diodes (OLED) with scattering layers, the plane wave propagation model of OLEDs was revised and its ability to rigorously simulate all electromagnetic properties of devices was demonstrated. The important exact relationships for the power flux and losses in OLED layers were given. The applicability of this model was confirmed by a comparison with the experimental data obtained by the measurement of the optical properties of fabricated green OLEDs.

Finally both methods for the light scattering calculation and OLEDs simulation were joined, and a scattering layer was demonstrated to increase an OLED external efficiency by several percent.

Introduction

Light-emittimg devices passed a long and exciting way from Edison's carbon glowers to modern diodes, and their development still face scientists with numerous interdisciplinary problems. Primary research motivators inlude the necessity of decreasing the cost of ligth while incresing its efficiency and flexibility. Lightning devices also suffer from modern trends of miniaturization, and it is hard colaborative work of scientists, engineers, and technicians that gives birth to the state-of-the art light-emitting diodes (LED).

Among the diversity of species, organic light emitting diodes (OLED) emerged after the 1987-th breakthrough of Tang and VanSlyke [START_REF] Tang | Oragnic electroluminescent diodes[END_REF] possess quite facinating potential. Their primary application sphere embraces large-area, possibly semi-transparent, thin ligtning panels and flexible displays as well as possible futuristic lightning decorations. Great effort was spent to bring them to the current almost ready-to-use state, and some further steps are required. This thesis presents an attempt to get ahead.

A history of OLEDs is a path from fractions to dozens of lm/W efficiency, and from several to many thousands of hours lifetime. Now the advance needs a solid foundation of sophisticated computer simulation. This work demonstrates results on rigorous numerical methods development for the light scattering and diffraction calculation in planar structures being representitative models for OLED optical properties study. Owing the powerful tool, it was applied to analyze a promising way of the OLED external efficiency increase by use of scattering layers. And, hopefully, this tool will also find its applications in other challenging problems.

The thesis is organized as follows. First chapter briefly describes a background of the treated problems including both review of concurrent numerical methods and description of problem area in the scope of organic light-emitting diode physics. Second chapter provides the descriprion of a new numerical method development key points as well as benchmarking results. Passage to the applications in OLED optics is made in the third chapter where an increase of OLED efficiency due to a scattering layer is demonstrated through the rigorous numerical simulation.

Present work is a result of collaboration between Laboratory of Nanooptics and Femtosecond Electronics in Moscow Institute of Physics and Technology, Dolgoprudny, Russia, and Laboratory Hubert Curien of University Jean Monnet of Saint-Etienne, France. I would like to thank my supervisors at both locations, Prof. Alexandre Tichshenko and Prof. Anatoly Gladun for invaluable assistance and help. Also, I wish to express my sencere gratitude to both collectives, in particular to Alexey Arsenin, who orchestrates the russian lab activity, and Prof. Olivier Parriaux in France. Finally, I cannot foget to mention my relatives, who I am grateful to for their kind empathy.

Chapter 1

Emission, propagation and scattering of light in planar structures

Electromagnetic wave propagation in homogeneous plane layered structures

Subject of this thesis is development and applications of models describing the light scattering and diffraction calculation in non-homogeneous planar dielectric structures. This work bases on a new numerical method which exibits less numerical complexity and computer memory requirements than known alternative approaches for the solution of a certain class of problems. Being developed the method is applied to a problem of rigorous simulation of organic ligh-emitting diodes (OLED) with scattering layers. Problem of the light diffraction and scattering in planar structures under consideration can be naturally divided into several sub-problems including the light propagation calculation in planar homogeneous structures, electroluminescent sources modeling, and light diffraction and scattering caluclation in spatially inhomogeneous layers. Each sub-problem was studied previously to some extent. Some known methods were taken as a basis for the current study. These methods together with those being close to the developed ones are briefly described in this chapter.

Generally a planar structure under consideration can be thought of a finite set of N L adjacent plane layers of different materials. These layers can be either homogeneous or inhomogeneous. Introduce Cartisian coordinates with axis Z being perpendicular to layers' plain, and designate the coordinates of plane interfaces between layers as z 0 . . . z N L (Fig. 1.1). Half-infinite media bounding the structure from above and from below with respect to axis Z positive direction will be further regarded as a cover and a substrate respectively. Denote coordinates of their boundaries with the structure as z L = z 0 и z U = z N L . Additionally, layer thicknesses are found as Emission, propagation and scattering of electromagnetic waves are considered here within the framework of classical electrodynamics, and, hence, are described by the Maxwell's equations:

h k = z k+1 -z k with k ∈ Z + : 0 ≤ k < N L . z y z 0 z 1 z N L z N L -1 h 0 h N L -1 ε 0 ε s ε c ε N L -1 a + c a - c a +
∇ × E = - ∂B ∂t , (1.1)
∇ × H = J + ∂D ∂t , (1.2)
∇ • D = ρ e , (1.3)
∇ • B = 0. (1.4)
In what follows all the fields and sources are considered to be decomposable into a set of monochromatic time harmonics with frequency ω, and each harmonic has an exponential factor of exp(-iωt). Assuming an external charge density to be 0 rewrite Eqs. (1.1)-(1.4) as a system of two equations for unknown time harmonic amplitudes:

∇ × E = iωB, (1.5)
∇ × H = J -iωD.

(1.6)

Eqs. (1.5) and (1.6) should be also supplemented with material relations between fields and inductions as well as electromagnetic boundary conditions [START_REF] Born | Principles of optics[END_REF]. We introduce these boundary conditions in a standard form of a matrix-vector product for the electric field, and in form of scalar product for the magnetic field:

D = εE, B = µ 0 H. (1.7)
With a view of writing out the explicit form of the first relation in (1.7) consider functions

ε k (r), k ∈ Z + : k = 0 . . . N L-1 , z k < z < z k+1
, describing the dielectric permittivty of the all layers in the planar structure (see Fig. 1.1). Here r = (x, y, z) is a vector in the intro-duced reference frame. Additional constants ε s and ε c correspond to substrate and cover spatially homogeneous permittivities. Composite layers consisting of different materials are described by discontinuous functions ε k (r). Therefore, the boundary conditions of the E and H vectors tangential components', and D and B normal components' continuity must be formulated at both set of interfaces z = z k , k = 0 . . . N L and surfaces of functions'

ε k (r) discontinuities.
Maxwell's equations (1.5), (1.6) contain external source density J. Here we single out two types of sources that will be considered further. The first one includes infinitely distant sources which radiation comes to a structure in form of plane waves. These waves are characterized by two wavevector projections k x , k y on axes X and Y , and their electric field amplitude writes:

E inc (r) = E inc exp(ik x x + ik y y ± ik z z), (1.8)
where

k z = ω 2 µ 0 ε -k 2 x -k 2 y .
(1.9)

The wavevector of a plane wave will be further regarded as k ± = (k x , k y , ±k z ). Here the sigh "+" distinguishes waves propagating in positive direction with repect to axis Z, whereas the sign "-" in negative. The second type of sources comprises dipole sources placed inside a structure or in the near-field region nearby it. A classical point monochromatic dipole placed at point r 0 = (x 0 , y 0 , z 0 ) is described in terms of the dipole moment density p(r, t) = p 0 δ(x -x 0 )δ(y -y 0 )δ(z -z 0 ) exp(-iωt).

(1.10)

Decomposing the electric and magnetic fields into sets of plane waves

f (r, t) = 1 (2πk 0 ) 2 ∞ -∞ ∞ -∞ f (k x , k y , t)dk x dk y , (1.11)
with vector f standing for both fields E and H, and factor k 0 = ω √ µ 0 keeps the dimensions, we come to the following formulas for the fields emitted by point dipole (1.10) (e.g., [START_REF] Polerecky | Theory of the radiation of dipoles placed within a multilayer system[END_REF])

E p (k x , k y , t) = k 2 0 2iε m k z k ± × k ± × p 0 exp ik ± r -iωt , H p (k x , k y , t) = iωk 2 2k z k ± × p 0 exp ik ± r -iωt .
(1.12)

Here ε m is the dielectric permittivity of a homogeneous isotropic medium containing the dipole source. In a series of works studying the luminescent molecules placed in the vicinity of plane interfaces, two-dimensional Fourier transform (1.11) is replaced with the Fourier-Bessel transform of variable k ρ = k 2 x + k 2 y . The corresponding formulas for the dipole field components are found, e.g., in [START_REF] Kong | Electromagnetic fields due to dipole antennas over stratified anisotrpic media[END_REF][START_REF] Chance | Molecular fluorescence and energy transfer near interfaces[END_REF][START_REF] Novotny | Allowed and forbidden light in near-field optics. i. a single dipolar light source[END_REF].

It is convenient to study polarization properties of planar structures by decomposing the fields into two independent TE and TM polarizations. For each plane harmonic these polarizarions are defined with respect to the plane of incidence defined by axis Z and wavevector k. The elelctric field of the TE-polarized wave is perpendicular to this plane whereas for the TM polarization the perpendicular vector is H. The relation between field amplitudes and plane polarized harmonic amplitudes a e± и a h± are given in the Appendix A in general matrix-vector form which is used throughout the thesis.

In the simplest case all the layers of a planar structure are homogeneous and isotropic. For zero structure thickness the problem reduces to a plane interface separating two media which is described by the well-known Fresnel coefficients for plane wave reflection and transmission [START_REF] Born | Principles of optics[END_REF]. If a structure contains a few layers of nonzero thickness, its reflection and refraction coefficients can also be written in an explicit analytic form (for example, [START_REF] Wait | Electromagnetic waves in stratified media[END_REF]). In general case of a multilayer stack one can distinguish two approaches for the rigorous electromagnetic field calculation S-matrix [START_REF] Ko | Matrix method for tunneling in heterostructures: Resonant tunneling in multilayer systems[END_REF] and T-matrix [START_REF] Born | Principles of optics[END_REF][START_REF] Heavens | Optical Properties of thin films[END_REF] S-matrix is different and relates incoming and outgoing wave amplitudes:

a - s a + c = S 00 S 01 S 10 S 11 a + s a - c , (1.14)
and corresponts to the quntum mechanical scattering operator [START_REF] Amrein | Scattering theory in quanum mechanics[END_REF], which translates an intial state of a system to its final state. T-matrices represent a convenient tool for theoretical analysis of planar media (e.g., [START_REF] Pelster | Propagation of plane waves and of waveguide modes in quasiperiodic dielectric heterostructures[END_REF][START_REF] Katsidis | General transfer-matrix method for optical multilayer systems with coherent, partially coherent, and incoherent interference[END_REF]) since their multiplication rule coincides with the standard matrix multiplication. However, in a numerical implementation T-matrices are unstabe allowing for the exponential error accumation in calculation of evanescent waves propagation [START_REF] Preist | Periodic multilayer gratings of arbitrary shape[END_REF]. On the other hand S-matrix-based numerical methods are stable, while their multiplication for a multilater structures is not so trivial. With a view of developing numerical methods we will deal only with S-matrices. Analytical expressions for S-matrix components have the simplest form for a plane interface between two media and for a plane homogeneous layer. An explicit form of the S-matrix for a plane interface separating two homogeneous media with dielectric permittivities ε L and ε U (Fig. 1.2a) writes via Fresnel coefficients [START_REF] Born | Principles of optics[END_REF]:

S T E I =     k L z -k U z k L z + k U z 2k U z k L z + k U z 2k L z k L z + k U z k U z -k L z k L z + k U z     (1.15)
for the TE polarization and

S T M I =     ε U k L z -ε L k U z ε U k L z + ε L k U z 2ε L k U z ε U k L z + ε L k U z 2ε U k L z ε U k L z + ε L k U z ε L k U z -ε U k L z ε U k L z + ε L k U z     (1.16)
for the TM polarization. S-matrix of a plane homogeneous layer of thickness h does not depend on the polarization state (Fig. 1.2b):

S L = 0 exp(ik z h) exp(ik z h) 0 . (1.17)
In many problems it is more convenient to deal with a set of plane harmonics rather than with a separate one. In this case S-matrix size equals to If S L and S U are written for a set of harmonics, equations (1.18) are matrix equations with the division operation being the multiplication by an inverse matrix. Notice also that operations sequence in (1.18) remains correct for S-matrix components S mn being either scalars or matrces. Equations (1.15)-(1.18) allow to calculate wave propagation, reflection and refraction in a planar structure for the light emitted by infitely distant sources. To simulate the emission of sources placed inside a structure decompose the field emitted by the point dipole (1.12) intoa set of TE and TM polarized waves [START_REF] Polerecky | Theory of the radiation of dipoles placed within a multilayer system[END_REF]:

2N O × 2N O instead of
a e± d = iω 2 µ 0 k 2 0 8π 2 γk z (k y p 0x -k x p 0y ) , a h± d = iω 2 k 2 0 8π 2 γk z ∓k x k z p 0x ∓ k y k z p 0y + γ 2 p 0z , (1.19)
where γ = k 2 x + k 2 y . Analogous formulas were combined with the T-matrix method to simulate the fluorescent molecule lifetine near plane interfaces in [START_REF] Kong | Electromagnetic fields due to dipole antennas over stratified anisotrpic media[END_REF][START_REF] Lukosz | Light emission by magnetic and electric dipoles close to a plane interface. i. total radiated power[END_REF][START_REF] Lukosz | Light emission by magnetic and electric dipoles close to a plane interface. ii. radiation patterns of perpendicular oriented dipoles[END_REF][START_REF] Lukosz | Theory of optical-environment-dependent spontaneous emission rates for emitters in thin layers[END_REF][START_REF] Crawford | Radiation from oscillatind dipoles embedded in a layered system[END_REF][START_REF] Benisty | Method of source terms for dipole emission modifcation in modes of arbitrary planar structures[END_REF][START_REF] Polerecky | Theory of the radiation of dipoles placed within a multilayer system[END_REF]. Similar approaches were used for the field decomposition into a set of cylindrical harmonics in [START_REF] Chance | Molecular fluorescence and energy transfer near interfaces[END_REF][START_REF] Wasey | Effects of dipole orientation and birefringence on the optical emission from thin films[END_REF][START_REF] Danz | Dipole lifetime in stratifed media[END_REF][START_REF] Ruppin | Lifetime of an emitting dipole near various types of interfaces including magnetic and negative refractive materials[END_REF]. Analytical methods based on Green's functions of a layered medium were developed in [START_REF] Yin | The theory of dyadic green's function and the radiation caracteristics of sources in stratified bi-isotropic media[END_REF][START_REF] Hartman | Green dyadic calculations for inhomogeneous optical media[END_REF][START_REF] Paulus | Accurate and effcient computation of the green's tensor for stratifed media[END_REF][START_REF] Hanson | Dyadic green's function for a multilayered planar medium a dyadic eigenfunction approach[END_REF][START_REF] Dogan | Closed-form representations of feld components of fluorescent emitters in layered media[END_REF], however, they appear to be interesting from the theoretical point of view, whereas in numerical computations they either allow obtaining only approximate results or reduce to T-matrix based calculations.

Diffraction and scattering in plane layers

The next sub-problem to be solved to develop a rigorous model of the light propagation and scattering in planar structures is calcuation of an inhomogeneous layer S-matrix. In the previous section S-matrices of homogeneous layers and plane interfaces we shown to have a quite simple form (1.15)-(1.17). However, for an inhomogeneous layer it seems to be impossible to obtain a general closed analytic form of the S-matrix. Thus, it has to be found numerically. Currently one can find several methods potentially capable to solve this problem. They can be classified into finite-difference, finite-element, integral equations, modal and hybrid methods.

Finite-difference (FD) methods form a broad class of methods capable to solve various differential equations. The core idea is the replacement of derivatives in some differential equation Lf = u by finite differences on a mesh G defined in a region x ∈ D where one searches for a solution f (x) ∈ S. The initial equation is replaced then by a finitedifference one L g f g (x) = u g which solution belongs to a mesh functions space S g . Principal properties of a finite-difference scheme include convergence, approximation and stability [START_REF] Bakhvalov | Numerical methods (in russian)[END_REF]. Convergence implies the decrease of the difference between the mesh solution and the exact solution projected to the mesh F g with the decrease of the mesh step τ proportionally to an integer power (convergence rate) of this step:

f g -F g ≤ const • τ k
Approximation shows the precision of the mesh equation with the exact solution being substituted in it (this is called residual):

L g F g -u g ≤ const • τ k
Stability means that small perturbations in the initial data lead to small changes in solution uniformly over the mesh step:

z ∈ S g , L g ϕ g = 0, L g ψ g = z ⇒ ϕ g -ψ g ≤ const • z
One of the most popular finite-difference schemes was proposed in [START_REF] Yee | Numerical solution of inital boundary value problems involving maxwell's equations in isotropic media[END_REF] (the article was cited more than seven thousands times). This scheme is characterized by the second order approximation over the space and time. Computer programms based on it were applied, e.g., for the modeling of diffractive optical elements [START_REF] Wong | Rigorous three-dimensional time-domain finitedifference electromagnetic simulation for photolithographic applications[END_REF][START_REF] Guerrieri | Massively parallel algorithms for scattering in optical lithography[END_REF], photonic crystals [START_REF] Painter | Defect modes of a two-dimensional photonic crystal in an optically thin dielectric slab[END_REF][START_REF] Boyd | Sensitive disk resonator photonic biosensor[END_REF][START_REF] Noda | Spontaneous-emission control by photonic crystals and nanocavities[END_REF], light scattering in non-periodic scattering media [START_REF] Mishchenko | Scattering, absorption, and emission of light by small particles[END_REF][START_REF] Jurgens | Finite-differences time-domain modeling of curved surfaces[END_REF]. The mentioned scheme is based on a cubic mesh as are many others in the FD method. Obviously they do not suit well to problems with a complex geometry which needs an adaptive mesh generation to take into account a specific field distribution. Examples of non-uniform meshes for the FD method were also proposed but generally they substantially complicate the method's formulation [START_REF] Beilenhoff | Improved finite-difference formulation in frequency domain for three-dimensional scattering problems[END_REF]. In this sence the finite element method (FEM) much better describes problems with complex cuved interfaces.

For the problems of the Maxwell's equations solution, application of the FD method results in a sparse linear algebraic equation system. The complexity of standard iterative methods [START_REF] Tyrtyshnikov | Methods of numerical analysis[END_REF] is O(N 2 ) with N being the size of the equation system, however, the sparseness usually allows to reduce the complexity down to O(N 1+α ), α > 0 [START_REF] Kunz | Finite difference time domain method for electromagnetic[END_REF]. A substantial drawback of the method is the necessity of taking the dimensions of a computational domain several times larger than the dimensions of a scattering object. This results in a tremendous increase of the mesh node number and, consequently, to the increase of the required computer memory. Thus, application of the FD method to complex scatter-ing structures usually requires the use of computing clusters with singnificant amount of memory.

The finite difference method was widely used for the solution of Maxwell equation as well as numerous other differential equations systems due to its universality. As was mentioned above it faces certain difficulties for the problems where the boundary conditions should be defined on complex shaped surfaces. To avoid facing this problem and retaind the universality property one may choose the finite element method (in the electrodynamics it is often referred to as the method of moments MOM). The idea lying in the core of the FEM is to calulate a decomposition of the solution f (x) of a differential equation system over a complete set of orthogonal functions. Denoting an approximate solution as f (x) one can write such decomposition in form

f (x) = ψ 0 + n c n ψ n (x).
Coeficients c n are found in scope of either variational principle or minimal energy condition. The variational formulation implies the calculation of such function f ∈ Ŝ that ( f ′ , g ′ ) = (u, g) for any g ∈ Ŝ. The minimal energy condition requires the minization of the funtional

E(g) = 1 2 (g ′ , g ′ ) -(u, g): f ∈ S, E( f ) ≤ E(g) for any g ∈ Ŝ.
Both of these formulations yield linear system equations with sparce matrices. These systems are solved by iterative algorithms [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. Apart from the mesh generation problem, one of the main shortcomings of the method for complex 3D scattering structure analysis is the same as that of the FD method very high requirements to the computer memory and, consequently, quite time-consuming calculations. One may regard this problem to be a consequence of the universality.

Certain scattering structures can be also analyzed by means of the volume integral equation (VIE) and the surface integral equation methods [START_REF] Kahnert | Numerical methods in electromagnetic scattering theory[END_REF]. Consider a scattering volume V . The volume integral equation in the 3D coordinate space writes [START_REF] Martin | Electromagnetic scttering in polarizable backgrownds[END_REF]

E sca (r) = E inc (r) + V dr ′ k 2 m ∆ε(r ′ )G m (r, r ′ )E sca (r ′ ), (1.20)
where ∆ε is the difference in the dielectric permittivity of a scatterer and a surrounding medium ε m , k m = ω √ ε m µ 0 , and G m (r, r ′ ) is the free space tensor Green's function [START_REF] Morse | Methods of theoretical physics[END_REF].

Eq. (1.20) is reduced to an algebraic linear equation system by subdividing the volume V into a number of sub-volumes ∆V which dimensions are small comparable to the wavelength λ (a conventional estimation of the maximum sub-volume size is λ/20). A similar equation system arises in the discrete dipole approximation (DDA) [START_REF] Yurkin | The discrete dipole approximation: An overview and recent developments[END_REF]. These methods are promising in terms of the numerical complexity which is linear with respect to the number of spatial sub-volumes. However, the VIE and the DDA are generally restricted to simulation of single particles or small groups of particles with dimensions are comparable to the wavelength and are hardly applicable to high-aperture structures. Another integral method is the surface integral method (SI) which is also referred to as the null-field method (NFM) [START_REF] Kahnert | Numerical methods in electromagnetic scattering theory[END_REF]. It is based on the derivation of surface integral equations from the Maxwel's equations and further field decomposition into a set of spherical harmonics. The main application of this method is the T-matrix calculation of non-spherical bodies for further use in the T-matrix method [START_REF] Mishchenko | T-matrix theory of electromagnetic scattering by particles and its applications: A comprehensive reference database[END_REF] (note, that in the light scattering theory T-matrices are different from those given in (1.13) and substantially represent element S 11 of S-matrix (1.14)).

A general observation concerning numerical methods consists in higher preference of narrow specialized methods in comparison with widely applicable approaches such as FD and FEM. The most effective among all seem to be modal methods. They include, e.g., the Mie solution [START_REF] Mie | Beiträge zur optik trüber medien, speziell kolloidaler metallösungen[END_REF][START_REF] Boren | Absorption and scattering of light by small particles[END_REF] describing a plane wave scattering on a sphere, the Tmatrix method of the light scattering calculation [START_REF] Waterman | Symmetry, unitarity, and geometry in electromagnetic scattering[END_REF][START_REF] Peterson | T matrix for electromagnetic scattering from an arbitrary number of scatterers and representation of e(3) *[END_REF][START_REF] Borghese | Multiple electromagnetic scattering from a cluster of spheres. i. theory[END_REF][START_REF] Mishchenko | T-matrix theory of electromagnetic scattering by particles and its applications: A comprehensive reference database[END_REF], and the modal method of the diffraction calculation on gratings [START_REF] Petit | Electromagnetic theory of graings[END_REF][START_REF] Peng | Theory of periodic dielect waveguides[END_REF][START_REF] Botten | The dielectric lamellar diffraction grating[END_REF][START_REF] Botten | The finitely conducting lamellar diffraction grating[END_REF][START_REF] Botten | Highly conducting lamellar diffraction grating[END_REF][START_REF] Tayeb | On the numerical study of deep conducting lamellar diffraction gratings[END_REF]. In problems where it is possible to analytically represent the modal fields, modal methods demonstrate the precision, speed and convergence far better than all other methods (for example, the light scattering calculation by a group of spherical particles or the light diffraction calculation on a rectangular grating). The other side of the matter is that an analytical representation highly narrows the range of modal methods direct applicability.

A way to improve modal methods' capabilities is to use a transformation to the Fourier space. This is done in the Fourier-modal method (FMM) [START_REF] Knop | Rigorous diffraction theory for transmission phase gratings with deep rectangular grooves[END_REF] also referred to as the rigorous coupled-wave approach (RCWA) [START_REF] Moharam | Rigorous coupled-wave analysis of planargrating diffraction[END_REF], and in the differential method. The FMM is widely used for the light diffraction calculation on plane gratings and diffractive optical elements. These methods are descibed below in some more detail than others since the method developed in this work also operates in the Fourier space. Formely the FMM was developed for the light diffraction caluclation on gratings so it will be described from this point of view.

Consider a plane periodically structured plane layer (Fig. 1.3 shows an example of a 2D sinusoidal intefrace separating two media within a plane layer). Such structuring is described as periodic change of the dielectric (and, perhaps, magnetic) permittivity along one or two noncollinear directions in plane XY . Let Λ 1,2 be the grating periods and k1,2 are unit vectors in the directions of periodicity. Then, material constants of the layer are written as spatial coordinate functions:

ε(r) = ε(r + m 1 Λ 1 k1 + m 2 Λ 2 k2 ), µ(r) = µ(r + m 1 Λ 1 k1 + m 2 Λ 2 k2 ), (1.21)
with integers m 1,2 , and

K 1,2 = k1,2 K 1,2 , K 1,2 = 2π/Λ 1,2 .
Decompose the electromagnetic field in the grating layer into a set of spatial harmonics:

f (r) = ∞ n 1 =-∞ ∞ n 2 =-∞ f n 1 n 2 (z) exp(in 1 K 1 ρ ρ ρ + in 2 K 2 ρ ρ ρ), (1.22)
where ρ ρ ρ = (x, y) is the radius-vector in plane XY , and f , stands for both the electric and the magnetic field. Indices n 1,2 enumerate diffraction orders. An inverse transform writes

f n 1 n 2 (z) = 1 Λ 1 Λ 2 Λ 1 0 Λ 2 0 exp(-in 1 K 1 ξ 1 -in 2 k 2 ξ 2 )f (ξ 1 , ξ 2 , z)dξ 1 dξ 2 .
(1.23)

Here ξ 1,2 are the coordinates in a frame which two axes Ξ 1,2 are defined by the reciprocal lattice vectors K 1,2 , and the third one coinsides with Z. Transformations (1.22) and (1.23) allow one to rewrite Maxwell's equations (1.5) and (1.6) in form of an infinite linear differential equation system with functions depending on the variable z:

dE x n 1 n 2 dz = iωµ 0 H y n 1 n 2 - k x n 1 ,n 2 ω m 1 ,m 2 ε (n 1 -m 1 )(n 2 -m 2 ) k x m 1 m 2 H y m 1 m 2 -k y m 1 m 2 H x m 1 m 2 , dE y n 1 n 2 dz = -iωµ 0 H x n 1 n 2 - k y n 1 ,n 2 ω m 1 ,m 2 ε (n 1 -m 1 )(n 2 -m 2 ) k x m 1 m 2 H y m 1 m 2 -k y m 1 m 2 H x m 1 m 2 , dH x n 1 n 2 dz = -iω m 1 ,m 2 ε (n 1 -m 1 )(n 2 -m 2 ) E y m 1 m 2 + k x n 1 ,n 2 ωµ 0 k x n 1 n 2 E y n 1 n 2 -k y n 1 n 2 E x n 1 n 2 , dH y n 1 n 2 dz = iω m 1 ,m 2 ε (n 1 -m 1 )(n 2 -m 2 ) E y m 1 m 2 + k x n 1 ,n 2 ωµ 0 k x n 1 n 2 E y n 1 n 2 -k y n 1 n 2 E x n 1 n 2 .
(1.24)

This system contains the Fourier-images of the dielectric permittivity ε n 1 n 2 defined in accordance with (1.23), as well as the Fourier images of the inverce permittivity:

ε n 1 n 2 (z) = 1 Λ 1 Λ 2 Λ 1 0 Λ 2 0 1 ε(ξ 1 , ξ 2 , z) exp(-in 1 K 1 ξ 1 -in 2 k 2 ξ 2 )dξ 1 dξ 2 . (1.25)
Wavevector projections for different diffraction orders k α n 1 n 2 , α = x, y, are defined as

k α n 1 n 2 = k inc α + n 1 K 1α + n 2 K 2α , α = x, y, (1.26)
where k inc α are the in-plane projections of the incident plane wave wavevector (1.8). Next, we rewrite the differential equation system (1.23) in matrix form

F ′ (z) = HF, (1.27)
where vector F = (E x , E y , H x , H y ) contains all the field harmonic amplitudes. The foundation of the FMM was laid in works [START_REF] Burckhardt | Diffraction of a plane wave at a sinusoidally stratified dielectric grating[END_REF][START_REF] Knop | Rigorous diffraction theory for transmission phase gratings with deep rectangular grooves[END_REF][START_REF] Moharam | Rigorous coupled-wave analysis of planargrating diffraction[END_REF]. This method consists in search of eigen solutions of Eq. (1.27), or, in other words, search of modes in the reciprocal Fourier space. First, a grating layer is divided into a number of thin slices (slicing approximation), and in each slice the dependence of the ε from coordinate z is neglected. Then the equation

HF = ±βF (1.28)
is used to find propagation constants of the Fourier modes. For numerical calculations the size of matrix H is made finite by cutting the infinite Fourier series. Denote the corresponding maximum diffraction order numbers as max |n 1,2 | = N O1,2 . Then Eq. (1.28) becomes a matrix equation with matrix H N O1 N O2 , and the problem is reduced to the algebraic eigenvalue problem. Numerical complexity of the last problem is

O(N 3 O1 N 3 O2 ) [37].
In the differential method one applies a finite-difference scheme to solve (1.27). There were proposed different implementations of this method, e.g., [START_REF] Petit | Electromagnetic theory of graings[END_REF][START_REF] Loewen | Diffraction gratings and applications[END_REF]. Currently the differential method is not as widely used as the FMM and have been developed in recent years mainly by authors of [START_REF] Nevière | Analysis of dielectric gratings of arbitrary profiles and thicknesses: comment[END_REF].

The above formulation of the Fourier-methods works well only for holographic gratings descibed by continuous functions (1.21). For corrugated one-dimensional gratings the corresponding method has a very poor convergence for TM waves diffraction. Authors of [START_REF] Granet | Efficient implementation of the coupled-wave method for metallic lamellar gratings in tm polarization[END_REF][START_REF] Lalanne | Highly improved convergence of the coupled-wave method for tm polarization[END_REF] demonstrated this problem to be caused by an incorrect passage from (1.27) to a finite system of linear equations. The essence of the problem is that the Fourier transform of the product of two functions having coincident points of discontinuity is incorrect. An attempt of mathematical justification of this fact was undertaken in [START_REF] Li | Use of fourier series in the analysis of discontinuous periodic structures[END_REF], however, one can think of a simpler explanation it is impossible to define a product of corresponding distributions [START_REF] Gushchin | Fourier modal method for relief gratings with oblique boundary conditions[END_REF][START_REF] Shcherbakov | Fast and memory-sparing exact electromagnetic analysis of arbitrary profile 2d periodic dielectric structures[END_REF][START_REF] Vladimirov | Generalized functions in mathematical physics[END_REF].

Appearence of works [START_REF] Granet | Efficient implementation of the coupled-wave method for metallic lamellar gratings in tm polarization[END_REF][START_REF] Lalanne | Highly improved convergence of the coupled-wave method for tm polarization[END_REF] stimulated both intense development and application of the Fourier-methods [START_REF] Granet | Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution[END_REF][START_REF] Popov | Grating theory: new equations in fourier space leading to fast converging results for tm polarization[END_REF], in particular, for the light difraction calculation on 2D gratings [START_REF] Lalanne | Improved formulation of the coupled-wave method for two-dimensional gratings[END_REF][START_REF] Popov | Maxwell equations in fourier space: fast-converging formulation for diffraction by arbitrary shaped, periodic, anisotropic media[END_REF][START_REF] Granet | Parametric formulation of the fourier modal method for crossed surface-relief gratings[END_REF][START_REF] David | Fast factorization rule and plane-wave expansion method for two-dimensional photonic crystals with arbitrary hole-shape[END_REF][START_REF] Schuster | Normal vector method for convergence improvement using the rcwa for crossed gratings[END_REF][START_REF] Götz | Normal vector method for the rcwa with automated vector field generation[END_REF]. Formulation of the FMM and the differential methods for 2D gratings requires an additional effort to treat the boundary condition in correct manner. Mathematical description of this problem will be given in the next chapter. There will be proposed a different approach from the one developed in [START_REF] David | Fast factorization rule and plane-wave expansion method for two-dimensional photonic crystals with arbitrary hole-shape[END_REF][START_REF] Schuster | Normal vector method for convergence improvement using the rcwa for crossed gratings[END_REF][START_REF] Götz | Normal vector method for the rcwa with automated vector field generation[END_REF].

Fourier-methods find their applications in the simulation and optimization of diffractive optical elements, photonic crystals (e.g., [START_REF] Hirayama | Rigorous electromagnetic analysis of diffractive cylindrical lenses[END_REF][START_REF] Lalanne | Design and fabrication of blazed binary diffractive elements with sampling periods smaller than the structural cutoff[END_REF][START_REF] Bezus | Calculation and modeling of diffractive structures for formation of 2d surface electromagnetic waves interference patterns[END_REF][START_REF] Turunen | Diffractive optics: Electromagnetic approach[END_REF]) which are examples of periodic and quasi-periodic structures. Recently there appeared several works where authors made attempts to adapt the FMM to the solution of non-periodic problems [START_REF] Edee | Complex coordinate implementation in the curvilinear coordinate method: application to plane-wave diffraction by nonperiodic rough surfaces[END_REF][START_REF] Pisarenco | Aperiodic fourier modal method in contrast-feld formulation for simulation of scattering from finite structures[END_REF][START_REF] Pisarenco | Modified s-matrix algorithm for the aperiodic fourier modal method in contrast-field formulation[END_REF]. These works describe calculaiton of the light diffraction on gratings containing the perfectly-matched layer (PML) [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Chew | A 3d perfectly matched medium from modified maxwell's equations with stretched coordinates[END_REF]. The PML allowed eliminating the re-scattering process on different grating periods and, hence, obtainig an approximate solution of a non-periodic problem. The formulation of the method developed in this thesis allows for the immediate incorporation of the PML, however, the simulation of scattering in this work will be carried out in a simpler manner.

Generalized source method

Now we proceed to description of a general theoretical method that was used in this thesis. This method was proposed in [START_REF] Tishchenko | A generalized source method for wave propagation[END_REF][START_REF] Tishchenko | Generalized source method: new possibilities for waveguide and grating problems[END_REF] where it was referred to as the generalized source method (GSM). GSM represents a rigorous procedure for calaultion of the light scattering and diffraction in inhomogeneous media. The method consists of two subsequent steps. First, one should choose a basis medium described by functions ε(r), µ(r) and a corresponding analytic solution of the Maxwell's equations for any source distribution. Second, this analytic solution shoud be written for a generalized source representing the difference between the initial and the basis media, thus, giving rize to a self-consistent equation. Now consider the described scheme in detail.

Let us start from the Maxwell's equations (1.5), (1.6), and rewrite the first one in a more general case allowing for the presence of inhomogeneous magnetic permittivity µ(r) and magnetic currents F(r):

∇ × E = iωµH + F.
(1.29)

Eqs. (1.29) and (1.5) give rize to the wave equations

∇ × 1 µ ∇ × E -ω 2 εE = iωεJ + ∇ × 1 µ F, (1.30)
∇ × 1 ε ∇ × H -ω 2 µH = -iωµF + ∇ × 1 ε J. (1.31)
Exact analytical solutions for (1.30) and (1.31) for any source distribution are known for a quite narrow range of problems. The GSM uses the power of these analytical solutions allowing one to develop numerical methods capable to solve wide classes of problems. To be specific, choose one of the exact solutions of system (1.30), (1.31) with certain boundary conditions and regard it to as the basis one. Write out the basis solution in form of a functional relationship that translates sources to unknown fields:

E = ℵ bE (J, F), H = ℵ bH (J, F). (1.32)
Decompose the permittivities describing an initial medium into a sum of the chosen basis permittivities and additional summands (which generally can be of any magnitudes): .33) This representation enables one to introduce generalized currents J gen and F gen generated by permittivity differences ∆ε(r) and ∆µ(r):

ε(r) = ε b (r) + ∆ε(r), µ(r) = µ b (r) + ∆µ(r). ( 1
J gen = -iω∆εE, (1.34)
F gen = iω∆µH, (1.35)
Then, Eqs. (1.32) can be rewritten in form

E = ℵ bE (J r + J gen , F r + F gen ) = E inc + ℵ bE (-iω∆εE, iω∆µH), H = ℵ bH (J r + J gen , F r + F gen ) = H inc + ℵ bH (-iω∆εE, iω∆µH). (1.36)
Here J r denotes real currents which are replaced by external fields E inc , H inc excited by them. Equation system (1.36) is a general form of implicit equations being the cornerstone of the GSM. Particular form of Eq. (1.36) depends on the basis solution. Note also that no restrictions are imposed to ∆ε(r), ∆µ(r).

The GSM can be thought of the following demonstrative procedure. Let an incident electromagnetic wave E inc , H inc be propagating in a region with scattering bodies. It excites the generalized sources which amplitudes are proportional to the functions ∆ε(r), ∆µ(r), and the modified field writes

E (1) = E inc + ℵ bE (-iω∆εE inc , iω∆µH inc ), H (1) = H inc + ℵ bH (-iω∆εE inc , iω∆µH inc ).
(1.37)

This field also interacts with the scattering structure and the subsequent modification is (1) , iω∆µH (1) ), (1) , iω∆µH (1) ).

E (1) = E inc + ℵ bE (-iω∆εE
H (1) = H inc + ℵ bH (-iω∆εE
(1.38)

And so on. Eq. (1.37) is, evidently, the Born approximation [START_REF] Landau | Quantum mechanics: non-relativistic theory[END_REF]. By continuing with the described iterations up to infinity we run into the Newmann series [START_REF] Tsang | Scattering of electromagnetic waves. Advanced topics[END_REF]. However, this treatment is good only for undertanding the method since in numerical calculations the Newmann series often diverge. For simulation of the light scattering and diffraction on high-constrast and high-aperture objects more sophisticated numerical methods should be used. The VIE method [START_REF] Martin | Electromagnetic scttering in polarizable backgrownds[END_REF] described above is an example of the GSM implementation with the basis medium being an isotropic homogeneous space and the basis solution being the tensor Green's function of a free space [START_REF] Akhmedzhanov | Modeling of the light scattering by nanoparticles of complex shape using the generalized source method[END_REF]. In this thesis the GSM is applied for the development of the method for the light scattering and diffraction calculation in plane micro-and nanostructured plane layers. An attempt of developing a similar method was undertaken in [START_REF] Magath | Fast iterative, coupled-integral-equation technique for inhomogeneous profiled and periodic slabs[END_REF][START_REF] Magath | Coupled integral equations for diffraction by profiled, anisotropic, periodic structures[END_REF], however authors succeded only in calculation of the TE wave diffraction on 1D gratings. A better result was obtained in [START_REF] Van Beurden | Fast convergence with spectral volume integral equation for crossed block-shaped gratings with improved material interface conditions[END_REF] for the rectangular crossed gratings calculation on basis of a method similar to the one proposed here. This article has appered recently and almost simultaneously with the article [START_REF] Shcherbakov | Fast and memory-sparing exact electromagnetic analysis of arbitrary profile 2d periodic dielectric structures[END_REF] describing the method proposed here.

Organic light emitting diodes with scattering layers

In this section we discuss the problem of rigorous simulation of organic light emittins diodes (OLED) with scattering layers, which is solved in the third chapter the thesis.

Electroluminescent properties of organic materials have been studied since about 50 th of the previous century. The first LED made of organic materials was created in 1989 [START_REF] Tang | Oragnic electroluminescent diodes[END_REF]. This breakthrough gave rise to a new field in science and technology which have been intensively developed so far. Currenlty the market offers a set of small-size OLED de-vices including displays for cellular phones and household appliances as well as decorative lightning elements. Leading companies announce the mass production of high-diagonal TV-sets and high-quality lightning devices for the next several years. Herewith manufacturers still face a range of problems including the prolongation of the OLED lifetime, search of new functional materials and increase of the device efficiency. The last problem is particularly important for lightning applications where OLEDs face quite a high competition with inorganic LEDs.

An OLED conventionally represents a multilayered structure showed in Fig. 1.4. A typical structure contains metal or transparent electrodes, and organic electroluminescent layer. Also there may be included additional organic layers allowing to tune the electronhole transport in a device and the device color. The OLED efficiency is defined as a ratio between the number of emitted photons and the number of electrons passed through a device. An alternative definition is the raio between the emitted light power and the electric power consumed by a diode. The power losses are divided into two essentially different channels: power loss due to the nonradiative exciton recombination with so called internal efficiency coefficient η in , and power loss due to optical trapping in the diode multilayer structure (including a substrate) described by the external efficiency coefficient η out . The net efficiency then writes

η ext = η in η out .
(1.39)

Currently the internal efficiency can be made close to 100% due to the use of phospho- rescent materials [START_REF] Adachi | Nearly 100% internal phosphorescence efficiency in an organic light emitting device[END_REF][START_REF] So | Recent progress in solution processable organic light emitting devices[END_REF][START_REF] Yersin | Highly efficient OLEDs with phosphorescent materials[END_REF]. Thus, the main effort regarding the solution of the OLED efficiency enhancement problem is directed to the improvement of η out . Optical losses can also be divided into several channels [START_REF] Nowy | Light extraction and optical loss mechanisms in organic light-emitting diodes[END_REF]. First, the power is absorbed in OLED layers, mainly in a metal cathode (up to 50% of the electromagnetic power). Second a substantial amount of power can be guided in the waveguide modes (∼10%). Additionally the power loss occures at the substrate-air interface due to the total internal reflection (up to ∼10%).

To lower the losses due to the waveguide modes excitation and the total internal reflection there were proposed several methods [START_REF] Meerholz | Outsmarting waveguide losses in thin-film lightemitting diodes[END_REF][START_REF] Nowy | Light extraction and optical loss mechanisms in organic light-emitting diodes[END_REF][START_REF] Riedel | Methods for increasing the efficiency of organic light emitting diodes[END_REF][START_REF] Wei | Emission characteristics of organic light-emitting diodes and organic thin-films with planar and corrugated structures[END_REF] (Fig. 1.5) including the use of photonic crystals and diffraction gratings [START_REF] Boroditsky | Light extraction from optically pumped light-emitting diode by thin-slab photonic crystals[END_REF][START_REF] Lupton | Bragg scattering from periodically microstructured light emitting diodes[END_REF][START_REF] Revelli | Waveguide analysis of organic lightemitting diodes fabricated on surfaces with wavelength-scale periodic gratings[END_REF][START_REF] Kim | Nanohole-templated organic light-emitting diodes fabricated using laser-interfering lithography: moth-eye lighting[END_REF][START_REF] Ziebarth | A theoretical and experimental investigation of light extraction from polymer light-emitting diodes[END_REF][START_REF] Vandersteegen | Employing a 2d surface grating to improve light out coupling of a substrate emitting organic led[END_REF][START_REF] Geyer | Large-scale patterning of indium tin oxide electrodes for guided mode extraction from organic light-emitting diodes[END_REF][START_REF] Tutt | Distribution of radiation from organic light-emitting diode structures with wavelength-scale gratings as a function of azimuth and polar angles[END_REF][START_REF] Hauss | Periodic nanostructures fabricated by laser interference lithography for guided mode extraction in oleds[END_REF][START_REF] Cho | Solution-processed photonic crystals to enhance the light outcoupling efficiency of organic light-emitting diodes[END_REF], microlences [START_REF] Gu | Highexternal-quantum-efficiency organic light-emitting devices[END_REF][START_REF] Madigan | Improvement of output coupling efficiency of organic light-emitting diodes by backside substrate modification[END_REF][START_REF] Yamasaki | Organic light-emitting device with an ordered monolayer of silica microspheres as a scattering medium[END_REF][START_REF] Möller | Improved light out-coupling in organic light emitting diodes employing ordered microlens arrays[END_REF][START_REF] Peng | Coupling efficiency enhancement in organic light-emitting devices using microlens array theory and experiment[END_REF][START_REF] Lee | Efficiency improvement and image quality of organic light-emitting display by attaching cylindrical microlens arrays[END_REF][START_REF] Lin | Luminance and image quality analysis of an organic electroluminescent panel with a patterned microlens array attachment[END_REF][START_REF] Pan | Fabrication and design of various dimensions of multi-step ashperical microlens arrays for oled package[END_REF], scattering layers [START_REF] Shiang | Application of radiative transport theory to light extraction from organic light emitting diodes[END_REF][START_REF] Liu | Microcavity top-emitting organic light-emitting devices integrated with diffusers for simultaneous enhancement of efficiencies and viewing characteristics[END_REF], aerogel layers [START_REF] Tsutsui | Doubling coupling-out efficiency in organic light-emitting devices using a thin silica aerogel layer[END_REF], microstructurization of a substrate [START_REF] Schnitzer | 30% external quantum efficiency from surface textured, thin-film light-emitting diodes[END_REF][START_REF] Windisch | Light-emitting diodes with 31outcoupling of lateral waveguide modes[END_REF][START_REF] Matterson | Increased efficiency and controlled light output from a microstructured light-emitting diode[END_REF][START_REF] Nakanishi | Improvement of the light extraction efficiency of top-emitting organic light-emitting diodes by a two-dimensional diffraction layer fabricated using self-assembled nanoparticles[END_REF][START_REF] Riedel | Enhancing outcoupling efficiency of indium-tin-oxide-free organic light-emitting diodes via nanostructured high index layers[END_REF][START_REF] Chen | Light extraction from organic light-emitting diodes for lighting applications by sand-blasting substrates[END_REF][START_REF] Kim | The effects of nanopillar surface texturing on the photoluminescence of polymer films[END_REF][START_REF] Ho | Efficiency enhancement of flexible organic light-emitting devices by using antireflection nanopillars[END_REF], and microresonator geometries [START_REF] Tsutsui | Sharply directed emission in organic electroluminescent diodes with an optical-microcavity structure[END_REF][START_REF] Grüner | Emission enhancement in single-layer conjugated polymer microcavities[END_REF][START_REF] Tokito | Microcavity organic light-emitting diodes for strongly directed pure red, green, and blue emissions[END_REF][START_REF] Peng | High-efficiency microcavity top-emitting organic light-emitting diodes using silver anode[END_REF][START_REF] Lee | Cavity effects on light extraction in organic light emitting devices[END_REF][START_REF] Lu | Design of microcavity organic light emitting diodes with optimized electrical and optical performance[END_REF]. The best results regarding the increase of the external efficiency were obtained with periodic wavelength-scale structures and microlences. However, a periodical structurization leads to a strong anisotropy of the emitted radiation which is quite undesirable for many applications, and the use of microlences requires a rather high-cost technology. Nanostructured scattering layers seem to be the most prospective from the point of view of the compromiss between the efficiency increase potential and the device cost. It is natural to place it between the transparent electrode and the substrate as shown in Fig. 1.4 to scatter waveguide modes and simultaneously avoid affecting the electron-hole transport. Tuning the parameters of a scattering layer and the whole OLED requires the ability of their optical properties simulation. The accuracy of such simulation should be at least about 1% since the net expected effect of the scattering layer application is about 1-10%. Furthermore, the model should account for the evanescent wave scattering which generally requires a rigorous solution of the Maxwell's equations. So, it is seen that the problem of OLED with sattering layer simulation is quite sophisticated.

Optical properties of OLEd with homogeneous layers containing electroluminescent sources were simulated with T-matrix method or analytical reflection and transmission coefficient calculation combined with the dipole representation of sources [START_REF] Neyts | Simulation of light emission from thin-film microcavities[END_REF][START_REF] Danz | Dipole lifetime in stratifed media[END_REF][START_REF] Lee | Electrical and optical simulation of organic light-emitting devices with fluorescent dopant in the emitting layer[END_REF][START_REF] Chen | Efficient and rigorous modeling of light emission in planar multilayer organic light-emitting diodes[END_REF][START_REF] Epstein | Optical emission from organic lightemitting diodes[END_REF][START_REF] Epstein | Electromagnetic radiation from organic light-emitting diodes[END_REF]. Other models include an approximate integration in the plane wavevectors space [START_REF] Kahen | Rigorous optical modeling of multilayer organic light-emitting diode devices[END_REF][START_REF] Chen | Electromagnetic modeling of organic light-emitting devices[END_REF][START_REF] Celebi | Simplified calculation of dipole energy transport in a multilayer stack using dyadic green's functions[END_REF]. Some works presented a waveguide analysis of OLED multilayer structures [START_REF] Lu | Optimization of external coupling and light emission in organic light-emitting devices: modeling and experiment[END_REF][START_REF] Revelli | Waveguide analysis of organic lightemitting diodes fabricated on surfaces with wavelength-scale periodic gratings[END_REF]. Examples of the external efficiency optimization for OLEDs with homogeneous layers are found in [START_REF] Cui | Optimization of light extraction from oleds[END_REF][START_REF] Peng | Coupling efficiency enhancement in organic light-emitting devices using microlens array theory and experiment[END_REF][START_REF] Nowy | Simulation based optimization of lightoutcoupling in organic light-emitting diodes[END_REF][START_REF] Mladenovski | Exceptionally efficient organic light emitting devices using high refractive index substrates[END_REF][START_REF] Krummacher | Efficiency analysis of organic light-emitting diodes based on optical simulation[END_REF][START_REF] Nowy | Understanding losses in OLEDs: optical device simulation and electrical characterization using impedance spectroscopy[END_REF]. OLED wiht gratings were simulated in [START_REF] Kikuta | Estimation method for the light extraction efficiency of light-emitting elements with a rigorous grating diffraction theory[END_REF][START_REF] Vandersteegen | Employing a 2d surface grating to improve light out coupling of a substrate emitting organic led[END_REF][START_REF] Khoshnegar | Diffraction analysis of extraction efficiency for photonic crystal based white light emitting diodes[END_REF]. Besides, there were attempts to simulate microstructured OLEDs with FEM and FD methods [START_REF] Lee | A high-extraction-efficiency nanopatterned organic light-emitting diode[END_REF][START_REF] Do | Enhanced light extraction effciency from organic light emitting diodes by insertion of a two-dimensional photonic crystal structure[END_REF][START_REF] Lee | Far-field radiation of photonic crystal organic light-emitting diode[END_REF][START_REF] Sun | Organic light emitting devices with enhanced outcoupling via microlenses fabricated by imprint lithography[END_REF][START_REF] Yan | Enhancement of light extraction efficiency in oled with twodimensional photonic crystal slabs[END_REF][START_REF] Xu | Enhancement of light extraction efficiency in oled with two-dimensional photonic crystal slabs[END_REF][START_REF] Yu | Nano-imprint fabrication and light extraction simulation of photonic crystals on oled[END_REF][START_REF] Jeon | Far-field radiation of photonic crystal organic lightemitting diode[END_REF][START_REF] Jeon | Ultraviolet nanoimprinted polymer nanostructure for organic light emitting diode application[END_REF][START_REF] Cho | Weak-microcavity organic light-emitting diodes with improved light out-coupling[END_REF][START_REF] Williams | Finite element simulations of excitonic solar cells and organic light emitting diodes[END_REF][START_REF] Pohl | Fast field solver for the simulation of large-area oleds[END_REF][START_REF] Cui | Optimization of light extraction from oleds[END_REF][START_REF] Altuna | Corrugated organic light emitting diodes for enhanced light extraction[END_REF]]. Here we do not discuss these works as the drawbacks of the corresponding methods were outlined above. Approximate models based on the radiation transfer theory were proposed in [START_REF] Shiang | Application of radiative transport theory to light extraction from organic light emitting diodes[END_REF][START_REF] Heikkila | Light extraction limits in textured gan-ingan light-emitting diodes: Radiative transfer analysis[END_REF], however, this approach does not meet the posed requirements.

Conclusions

The first chapter gives the nesessary information conserning the problems solved in the thesis. It contains a brief review of numerical methods potentially concurrent to those ones developed in this work and capable to rigorously solve the Maxwell's equations in plane inhomogeneous layers, description of the generalized source method being the basis of the teoretical developments presented below, and a short discussion of the questions related to the OLED with scattering layers simulation problem. One can conclude that, first, currently there is a strong need in fast numerical methods in the light scattering theory capable to deal with complex structures, second, that the Fourier-methods are a promising choise for the problem of the light scattering in planar structures, and, third, that the problem of rigorous simulation of OLED with nanostructured scattering layers was not solved previously with suffcicent and controlled accuracy.

Chapter 2

Ligth diffraction on 2D diffraction gratings

Introduction

This chapter describes the implementation of the generalized source method in the 2D reciprocal space. First, we will obtain analytical formulas for the S-matrix components of an infinitely thin inhomogeneous layer, and, second, an implicit equation describing the light diffraction on gratings will be formulated. The end of the chapter demonstrates the convergence analysis of the proposed fast method.

An interrelation between the inital nonperiodic problem of the light scattering in plane heterogeneous layers and the light diffraction on gratings calculation problem can be established from the following considerations. An approach developed in this work bases on the planar geometry of layers independently of a shape of scattering particles placed inside a layer. Thus, a plane wave representation is used here as a natural representation of such geometry. Mathematically this representation is expressed as the 2D Fourier transform of the fields and permittivities in the XY plane. In accordance to the convolution theorem the Fourier image of a two functions product is a convolution of their Fourier images:

F(f • g) = F(f ) * F(g).
(2.1)

Thus, the Fourier transform converts products of permittivities by fields in wave equations (1.30) and (1.31) into corresponding convolutions. Any numerical calculation requires a discretization, which in this case is the discretization of the reciprocal space. Then, one can notice, that the convolution is represented by a product of a Toeplitz matrix by a vector providing that the mesh in the reciprocal space is equidistant. This product in turn can be calculated by the FFT as was mentioned in the previous chapter. Therefore, a method that solves the light diffraction problem by means of matrix-vector product operations with only Toeplitz and diagonal matrices can be implemented with the linear comlexity with respect to the mesh node number. The charge of the speed is the equidistancy of a mesh in the reciprocal space, and the corresponding periodicity in the coordinate space.

Basis solution

The theoretical analysis in this thesis bases on the generalized source method described in section 1.3. Here we consider the first step of the GSM applied to the diffraction problem, namely, the derivation of the basis solution. Starting from the Maxwell's equations one can write the Helmgoltz's equations in a homogeneous isotropic medium with permittivities ε b and µ b allowing for both electric and magnetic sources:

∇ (∇E b ) -∆E b -ω 2 ε b µ b E b = iωµ b J + ∇ × F, (2.2)
∇ (∇H b ) -∆H b -ω 2 ε b µ b H b = -iωε b F + ∇ × J. (2.3)
Introduce vector A E , A H and scalar ϕ E , ϕ H potentials as

E b = -∇ϕ E + iωA E - 1 ε b ∇ × A H , (2.4)
H b = ∇ϕ H + iωA H + 1 µ b ∇ × A E . (2.5)
Being substituted into Maxwell's equations (1.5), (1.6) the vector potentials can be shown to satisfy the following equations

∆A E + k 2 b A E = -µ b J, (2.6
)

∆A H + k 2 b A H = ε b F, (2.7)
providing that the Lorentz gauge is used [START_REF] Landau | The classical theory of fields[END_REF]:

ϕ E = ∇A E iωε b µ b , (2.8
)

ϕ H = - ∇A H iωε b µ b . (2.9)
Due to the further Fourier transform discussed in the introduction to the current chapter, we write sources in form of plane currents:

J F = j(z) f (z) exp (ik x x + ik y y) . (2.10)
Solutions of (2.6) and (2.7) for the specified form of sources can be written as an integral over coordinate z [START_REF] Tishchenko | Generalized source method: new possibilities for waveguide and grating problems[END_REF] A E (r)

A H (r) =i exp (ik x x + ik y y) 2k z × ∞ -∞ µ b j(z) -ε b m(z) exp [ik z (z -z ′ ) ξ (z -z ′ )] dz ′ , (2.11)
where ξ denotes the difference of two Heaviside θ-functions:

ξ (z -z ′ ) = θ (z -z ′ ) -θ (z ′ -z) = 1, z -z ′ ≥ 0 -1, z -z ′ < 0 , (2.12)
and wavevector z-projection k z is defined by 1.9. The fields are found from (2.4) and (2.5) and represent a superposition of plane waves propagating upwards and downwards with respect to axis Z together with an additional term proportional to the source amplitude:

E α = exp (ik x x + ik y y) δ αz iωε b j z (z) + β=x,y,z z -∞ j β (z ′ ) Y + βα ωε b -f β (z ′ )X + βα exp [ik z (z -z ′ )] dz ′ + β=x,y,z z -∞ j β (z ′ ) Y - βα ωε b -f β (z ′ )X - βα exp [-ik z (z -z ′ )] dz ′    , (2.13)
H α = -exp (ik x x + ik y y) δ αz iωµ b f z (z) - β=x,y,z z -∞ f β (z ′ ) Y + βα ωµ b + j β (z ′ ) X + βα exp [ik z (z -z ′ )] dz ′ - β=x,y,z z -∞ f β (z ′ ) Y - βα ωµ b + j β (z ′ ) X - βα exp [-ik z (z -z ′ )] dz ′    . (2.14)
Here indices α и β stand for spatial coordinates x, y, and z, matrix elements write

Y ± αβ = k ± α k ± β -k 2 δ αβ 2k z , (2.15)
X ± αβ = e αγβ k ± γ 2k z , (2.16)
and δ αβ , e αβγ are Kronecker symbol and absolutely antisymmetric tensor respectively. Eq. (2.13) and (2.14) provide the required basis solution of the GSM. Note that these formulas can be also derived using the method of Green's functions [START_REF] Morse | Methods of theoretical physics[END_REF][START_REF] Felsen | Radiation and scattering of waves[END_REF]. Consider an important case when sources exist only in plane z = 0 and their amplitudes (2.10) are represented by δ-functions of coordinate z:

j β (z) = j β δ(z), f β (z) = f β δ(z).
(2.17)

Then the firt summand in Eqs. (2.13), (2.14) gives a singular disturbance in the XY plane.

To account for the polarization state of the electromagnetic radiation we introduce a standard tranfromation of the field amplitudes to the amplitudes of TE-and TM-polarized waves [START_REF] Born | Principles of optics[END_REF]. Corresponding notations are given in Appendix A. The relation between the TE-and TM-waves amplitudes and the source amplitudes is found by substitutuin of (2.17) into (2.13) and (2.14) and taking into account Eqs. (A.7), (A.8):

a e± reg (j, f ) = ωµ 0 k x 2γk z j y - ωµ 0 k y 2γk z j x ± k x 2γ f x ± k y 2γ f y - γ 2k z f z , (2.18)
a h± reg (j, f ) = ωε 0 k y 2γk z f x - ωε 0 k x 2γk z f y ± k x 2γ j x ± k y 2γ j y - γ 2k z j z . (2.19)
In the source plane z = 0 these expressions should be supplemented with singular summands To simplify the further analysis and possibly improve the numerical behaviour of the method we introduce a modified field by susbstracting singular terms in the source region:

a e± = a e± reg + a δe = a e± reg + δ(z) j z iωε 0 , (2.20)
a h± = a h± reg + a δh = a h± reg -δ(z) f z iωµ 0 . ( 2
Ẽx,y = E x,y Ẽz = E z - j z iωε b , (2.22)
Hx,y = H x,y

Hz = H z + f z iωµ b , (2.23)
so that the modified fields are regular everywhere and can be decomposed into amplitudes (2.18) and (2.19) only.

S-matrix based diffraction calculation

Eqs. (2.18) and (2.19) are witten for a single source plane harmonic (2.10). In general case one should account for all possible harmonics of fields and currents. For this purpose we introduce index n enumerating the Fourier orders. Accounting for the further truncation of the Fourier series denote the maximum harmonic numbers as N O1 и N O2 . Then a one-to-one correspondence between n and indices enumeraing the diffraction orders along each direction of periodicity n 1 , n 2 can be established in form

n = n 1 N O2 + n 2 , -N O1,2 < n 1,2 < N O1,
ãe± m = ωµ b k xm 2γk zm j ym - ωµ b k ym 2γ m k zm j xm ± k xm 2γ m f xm ± k ym 2γ m f ym - γ 2k zm f zm , (2.24
) (2.28)

ãh± m = ωε b k ym 2γ m k z f xm - ωε b k xm 2γ m k zm f ym ± k xm 2γ m j xm ± k ym 2γ m j ym - γ m 2k zm j zm . ( 2
By introducing notations

∆ ε,µ x,ymn = (ε, µ) x,y (ε, µ) b mn -δ mn , ∆ ε,µ zmn = δ mn - (ε, µ) z (ε, µ) b mn , (2.29)
one can write

S ee±± mn = i 2 k b 2 k xm γ m k zm ∆ ε ymn k xn γ n + k b 2 k ym γ m k zm ∆ ε xmn k yn γ n ± m ± n k xm γ m ∆ µ xmn k xn k zn γ n ± m ± n k ym γ m ∆ µ ymn k yn k zn γ n + γ m k zm ∆ µ zmn γ n , (2.30)
S he±± mn = i 2 ωµ b ± n k xm γ m k zm ∆ ε ymn k yn k zn γ n -± n k ym γ m k zm ∆ ε xmn k xn k zn γ n ± m k xm γ m ∆ µ xmn k yn γ n ∓ m k ym γ m ∆ µ ymn k xn γ n , (2.31)
S eh±± mn = i 2 ωε b ± n k ym γ m k zm ∆ µ xmn k xn k zn γ n ∓ n k xm γ m k zm ∆ µ ymn k yn k zn γ n ∓ m k xm γ m ∆ ε xmn k yn γ n ± m k ym γ m ∆ ε ymn k xn γ n , (2.32)
S hh±± mn = i 2 k b 2 k ym γ m k zm ∆ µ xmn k yn γ n + k b 2 k xm γ m k zm ∆ µ ymn k xn γ n ± m ± n k xm γ m ∆ ε xmn k xn k zn γ n ± m ± n k ym γ m ∆ ε ymn k yn k zn γ n + γ m k zm ∆ ε zmn γ n , (2.33)
where signs ''±'' и ''∓'' with index n correspond to incident field harmonics, and with index m correspond to diffracted field harmonics. Analogously, one can obtain S-matrix components of an infinitely thin layer of a corrugated grating. They appear to be rather bulky and are given in Appendix B. Eqs. (2.30)-(2.33) and their modifications (B.2)-(B.10) allow one to calculate the light diffraction on gratings with the S-matrix multiplication rule (1.18). Since the derived equations describe the diffraction on an infinitely thin plane layer, calculation of the diffraction on a thick layer requires slicing of a grating along axis Z into a finite number of sufficiently thin slices, and calculation of S-matrices for each slice. Thus, an algorithm for the S-matrix-based diffraction calculation is formulated as follows:

1. Slicing of a layer of thickness h with a grating into N S slices of thickness ∆h; Steps 2 and 3 of the given algorithm will be retained in the method based on a linear algebraic equation system solution, and will be discussed further.

As can be noticed, the last step is the most computationally complex part of the given algorithm. According to Eqs. (1.18) this step requires inversion of matrices of size O N S ). There was written a program based on the given algorithm for light diffraction calculation on 1D gratings in both collinear and noncollinear cases. Set of input parameters include: the wavelenght of an incident plane wave, angle of incidence (two angles in noncollinear case), grating period and depth, grating profile, parameters N O , N S , and permittivites of all materials. The output includes all S-matrix complex components.

N O × N O ,
Figs. 2.1a and 2.2a demonstrate convergence of the method with the increase of the slice number for a plane wave diffraction with λ = 0.6328µm and incidence angle 10 • on 1D rectangular and sinusoidal corrugated gratings. Parameters of both gratings were taken to be Λ = 1µm, h = 0.5µm and refractive index contrast 1.5. The same convergence rate was revealed also for other types of gratings sinusoidal index gratings and gratings consisting of infinite cylinders, in the range of periods from 100 nm to 10 µm and in the range of depths from 10 nm to 5 µm. One may notice that the slice numbers used for the diffraction calculation on rectangular gratings greatly exceeds values of N S for the sinusoidal one. This is due to the independence of the rectangular grating profiles from coordinate z, which enables the power-law multiplication of S-matrices instead of linear subsequent multiplications.

Figs. 2.1a and 2.2a reveal that the convergence has the power-law dependence from the inverse number of slices starting from a sufficienly large N S depending on N O . Thus, one may suppose that the solution (S-matrix components) can be represented as a polynomial of the variable 1/N S in some neighbourhood of the point 1/N S = 0:

a m (1/N S ) = a m (0) + ∞ k=1 a (k) M 1 N S k , (2.34)
where a m (0) denotes the exact solution. Then applying the Lagrange's polynomial interpolation scheme [START_REF] Tyrtyshnikov | Methods of numerical analysis[END_REF] of power p and calculating a free term of the polynomial The method of the light diffraction calculation developed in this section is interesting mainly due to its novelty. The given estimation of its accuracy shows that it does not provide any advantages in comparison with the FMM. Nevertheless, it can be modified so as to reduce the diffraction problem to an implicit linear equation system mentioned in section 1.3. In other words this means that instead of a sequentional calculation of slice S-matrices and their subsequent multiplication all the wave amplitudes in each slice will be calculated self-consistently at a time.

a m (0) ≈ p+1 k=1 a m (1/N Sk ) q =k [1 -N Sq /N Sk ] (2

Diffraction on index gratings

In this section the proposed S-matrix based method of the light diffraction calculation is modified so as to reduce its numerical complexity and replace the S-matrix multiplications by self-consistent linear algebraic equation system solution.

We start with rewriting Eqs. (2.27) и (2.27) in matrix-vector form

j (E) αm j (M ) αm = n -iωε b V E αβmn 0 0 iωµ b V H αβmn Ẽβn Hβn , (2.36)
where we introduced block-diagonal matrices

V E = diag{∆ ε x , ∆ ε y , ∆ ε z }, V H = diag{∆ µ x , ∆ µ y , ∆ µ z }.
(2.37)

The last two matrices are diagonal with respect to coordinate indices α, β, and their components are defined by (2.29). Then, substitution of Eq. (2.36) into basis solutions (2.13) and (2.14), and using formulas of Appendix A, gives an implicit integral equation relating amplitudes of incident and diffracted harmonics:

a m (z) = ∞ -∞ ∞ n=-∞ α,β          R ± m (z, z ′ ) P E αm P H αm × V E αβmn (z ′ ) 0 0 V H αβmn (z ′ ) × Q E βn Q H βn a inc n (z ′ )          dz ′ , (2.38)
where

R ± m (z, z ′ ) = ζ [± (z -z ′ )] exp [±ik zm (z -z ′ )] , (2.39)
ζ (z) =      1, z > 0 1 2 , z = 0 0, z < 0 , (2.40)
and matrices P E,H write

P E m = 1 2       -ωµ b k ym /(2γ m k zm ) ωµ b k xm /(2γ m k zm ) 0 -ωµ b k ym /(2γ m k zm ) ωµ b k xm /(2γ m k zm ) 0 k xm /γ m k ym /γ m -γ m /k zm -k xm /γ m -k ym /γ m -γ m /k zm       , (2.41)
P H m = 1 2       k xm /γ m k ym /γ m -γ m /k zm -k xm /γ m -k ym /γ m -γ m /k zm ωε b k ym /(2γ m k zm ) -ωε b k xm /(2γ m k zm ) 0 ωε b k ym /(2γ m k zm ) -ωε b k xm /(2γ m k zm ) 0       . (2.42)
Vector a(z) = a e+ mn a e- mn a h+ mn a h- mn T contains all amplitudes of TE-and TM-harmonics propagating in a plane with coordinate z.

To pass from Eq. (2.38) to a system of algebraic equations we use a standard approximation of the integral by a finite sum over slices introduced in the previous section. Consider as before a layer bounded by planes z = z (L) and z = z (U ) , z (U ) -z (L) = h, and divide it into N S slices of equal thickness ∆h = h/N S . Let wave amplitudes a e,h± mp to be defined at centers of slices with coordinates z p , p = 0, . . . , (N S -1),

z p = z (L) + (p + 1/2) ∆h.
(2.43)

Then, Eq. (2.38) reduces to a linear system of algebraic equations

a mp = N S -1 q=0 ∞ n=-∞ α,β R ± mpq P αm V αβmnq Q βn , (2.44)
where indices p and q enumerate slices, and matrix elements, that were introduced before, are replaced by more compact notations including both electric and magnetic sub-parts. Then we truncate infinite sums over diffraction orders by some maximum value N O = N Ox N Oy . This brings us a finite equation system which writes in matrix form as a = RPVQa inc = Aa inc .

(2.45)

It is naturally to define the amplitude vector of an incident field at the boundaries of a grating layer by constants a e,h- inc (z (U ) ) and a e,h+ inc (z (L) ). The amplitudes in each slice a inc to be substituted in (2.45) then write

a e,h+ inc (z q ) a e,h- inc (z q ) = a e,h+ inc (z (L) ) exp [ik z0 ∆h(q -1/2)] a e,h- inc (z (U ) ) exp [ik z0 ∆h(N S + 1/2 -q)]
.

(2.46)

Note that here a layer with a grating is supposed to be placed in a homogeneous medium with constant permittivity ε b . Generalization of the following results to a case of arbitrary substrate and cover will be given further. Now, in accordance with the GSM (1.36), unknown diffracted field amplitudes are found via the solution of a self-consistent problem represented by a system of linear algebraic equations:

a = (I -RPVQ) -1 a inc = (I -A) -1 a inc , (2.47)
where I is the identity matrix

I = I αβmnpq = δ αβ δ mn δ pq . (2.48)
To calculate the amplitudes of diffraction orders propagating outwards of the layer boundary a out , one should apply Eq. (2.45) one more time, so that a out = a inc + TPVQ(I -RPVQ) -1 a inc .

(2.49)

The introduced matrix T coherently transforms unknown amplitudes of the diffracted waves into corresponding amplitudes at layer boundaries. Its components are found analogously to (2.46) and write

T (U ) nq = exp [ik zn ∆h(N S + 1/2 -q)], T (L) nq = exp [ik zn ∆h(q -1/2)].
(2.50)

Finally, to calculate unknown amplitudes of the electric and magnetic fields components inside a grating layer we use (A.7), (A.8), (2.22) and (2.23):

E(z q ) =    I 0 0 0 I 0 0 0 I -∆ ε z    Q E (I -A) -1 a inc , (2.51)
H(z q ) =    I 0 0 0 I 0 0 0 I -∆ µ z    Q H (I -A) -1 a inc .
(2.52)

Outside the grating layer the modified field (2.22), (2.23) coincides with the real field and can be found as the product of vector a out by matrices (A.7) and (A.8).

Thus, the obtained equations (2.47), (2.49), (2.51) and (2.52) fully describe the problem of the light diffraction calculation on holographic gratings with continuous functions ε(r) and µ(r). Now we proceed to a numerical algorithm.

Numerical algorithm

To describe a numerical algorithm of solution of Eqs. (2.47), (2.49), (2.51) and (2.52) we recall the considerations given in the beginning of the current chapter regarding the fast matrix-vector multiplication possibility in the Fourier space. To be able to use this technique, solution should be found by an iterative method. Besides, the use of an interative procedures is necessary since for big values of N O и N S matrices of the derived linear systems cannot be inverted by the direct matrix inversion. Complexity of the di-rect inversion algorithm is O(N 3 O N 3 S ), which bounds the possibility of its use on modern personal computers by values N O ∼ 100 and N S ∼ 10.

There are several methods capable to solve equations like (2.47), (2.49) iteratively [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. These equations have complex non-symmetric dence matrices. The most widely used methods in the light scattering theory [START_REF] Smith | The biconjugate gradient method for electromagnetic scattering[END_REF][START_REF] Kahnert | Numerical methods in electromagnetic scattering theory[END_REF] for similar systems are the biconjugate gradient method (Bi-CG) [START_REF] Saad | Iterative methods for sparse linear systems[END_REF], its modifications [START_REF] Van Der | Bi-cgstab: A fast and smoothly converging variant of bicg for the solution of nonsymmetric linear systems[END_REF][START_REF] Sleijpen | Bicgstab(l) and other hybrid bi-cg methods[END_REF] and the generalized minimal residual method (GMRES) [START_REF] Saad | Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. Use of these methods to solve (2.47), (2.49) revealed that the GMRES appears to be the most reliable in terms of convergence. Thus, its use will be implied in all subsequent numerical examples.

Consider in detail the structure of matrix A appearing in Eqs. (2.47), (2.49) and (2.51). This square matrix can be thought to consist of N S × N S blocks each containing

(2N O -1) × (2N O -1) sub-blocks of size 4 × 4.
In accordance with (2.47) matrix A is the product of four matrices A = RPVQ. Matrix R contains exponential factors describing the propagation of plane harmonics between different slices, and matrix elements R npq depend only on the slice index difference (p -q). This matrix is Toeplitz with repect to these spatial indices, or, in other words R is block-Toeplitz. Analogously, matrix elements V nmp depend on the difference (n -m), and V is block-Toeplitz with respect to the Fourier indices, while it is diagonal with respect to p, q. Matrices P and Q are block-diagonal regarding to inner blocks of size 3 × 4 и 4 × 3. Toeplitz matrix can be expanded to circulant matrix, and a product of the latter by a vector can be found by the FFT.

To summarize, the proposed numerical algorithm bases on the GMRES whith matrixvector products being calculated by the FFT. Numerical complexity of the algorithm is

O [N S N O log (N S N O )].
For large N O and N S the complexity is linear relative to the product N S N O which is much better than both complexities of the GMRES with usual matrixvector multiplication O (N 2 S N 2 O ) and of the considered S-matrix algorithm O (N S N 3 O ). Thus, the numerical algorithm is formulated as follows:

• Calculation of the incident field harmonics amplitudes in each slice according to (2.46).

• Calculation of the Fourier-images of the dielectric and magnetic permittivities and pre-caluclation of the FFT from the obtained matrices.

• Pre-calculation of the FFT from matrix R.

• Solution of (2.47) by the GMRES with the FFT using the FFT-pre-calulated matrices.

• Calculation of diffraction orders amplitudes by Eq. (2.49) or/and amplitudes of the field projections by the Eq.(2.51).

Diffraction on corrugated gratings

In section 2.4 we obtained Eqs. (2.47), (2.49) and (2.51) for calculation of the light diffraction on holographic gratings described by continuous functions ε(r) and µ(r). Numerous applications require the study of composite structures with sharp interfaces separating different media. In this case the presented solution of the Maxwell's equations is incorrect regarding to the problem of Fourier-methods mentioned in Chapter 1. In spite of the fact that from the physical point of view one may expect a relative proximity of the solution of (2.47), (2.49) to an exact one for all types of gratings, the derivation of correct equations is essential for obtaining a good convergence and controling the accuracy of the results.

We start with analysis of corrugated gratings from the observation that the generalized currents are proportional to the electric D and magnetic B induction:

J gen = -iω(D -D b ), F gen = iω(B -B b ).
(2.53)

Further only the presence of the generalized electric sources will be considered with the magnetic permittivity of both scattering medium and basis medium being equal to the vacuum permittivity µ 0 . Problem solution with the megnetic sources is absolutely analogous to the following analysis for electric sources. Thereby matrix indices introduced in Appendix A distinguishing electric and magnetic fields will be omitted since only the electric field will participate in the following derivations. As stated above, the correct treatment of generalized sources requires exclusion of discontinuous function products with coincidenting points of discontinuity. For this purpose we start with boundary conditions for the normal and the tangential components of the electric field at interfaces. The electric field tangential component E is continuous at interfaces, and taking the Fourier image in (2.53) is correct:

(D ) n = N O m=-N O [ε] nm (E ) m .
The normal component E ⊥ is discontinuous together with function ε(r) so that they must be separated at different parts of the material relation (1/ε)D = E. This results in

N O m=-N O [1/ε] nm (D ⊥ ) m = (E ⊥ ) n .

The last relation brings matrix-vector formula

D ⊥ = [1/ε] -1 E ⊥ to be used.
The stated relations together with (2.53) bring the normal and tangential components of the generalized current

(j ) n = -iωε b N O m=-N O ε ε b nm -I nm (E ) m , (j ⊥ ) n = -iωε b N O m=-N O ε b ε ⊥ -1 nm -I nm (E ⊥ ) m .
(2.54)

To apply these equations introduce a local coordinate system n, ψ and ϕ at a grating corrugation surface. Axis n coinsides with the normal direction to this surface. The other two axes are defined by angle ψ between normal direction and axis Z of the initial Cartesian coordiane system, and angle ϕ between the normal projection on plane XY and axis X. Axes ψ and ϕ ϕ ϕ lie in the tangential plane to a grating surface with axis ϕ ϕ ϕ lying in plane XY as shown in Fig. 2.3. For any vector b in the initial frame XY Z, its components in the local frame are found from the following transform:

   b n b ψ b ϕ    =    cos ϕ sin θ sin ϕ sin θ cos θ cos ϕ cos θ sin ϕ cos θ -sin θ -sin ϕ cos ϕ 0       b x b y b z    . (2.55)
An inverce transform writes via the transposed matrix as

   b x b y b z    =    cos ϕ sin θ cos ϕ cos θ -sin ϕ sin ϕ sin θ sin ϕ cos θ cos ϕ cos θ -sin θ 0       b n b ψ b ϕ    . (2.56)
Note, that the Jacobian of this transform equals to 1.

Next, suppose that the trigonometric functions of angles ψ and ϕ are smooth functions of coordinates x, y, z except, maybe of a finite number of points. This assumption is not strong since these finctions are initially defined on a set of curves in plane of each slice and they can be extended to the entire slice planes with periodic smooth functions with narrow spectra. The main condition to be preserved here is the requirement of main trigonometric identities validity on the whole slice plane:

sin 2 ψ(r) + cos 2 ψ(r) = 1, sin 2 ϕ(r) + cos 2 ϕ(r) = 1.
(2.57)

Approaches to the extrapolation of the mentioned trigonometric functions will be discussed below while considering particular grating examples. Analogous treatment of normal directions can be found in [START_REF] Gushchin | Fourier modal method for relief gratings with oblique boundary conditions[END_REF] relative to the FMM and in [START_REF] Popov | Grating theory: new equations in fourier space leading to fast converging results for tm polarization[END_REF] relative to the differential method. An alternative approach given in [START_REF] Schuster | Normal vector method for convergence improvement using the rcwa for crossed gratings[END_REF] in detail is to generate the normal vector field on a grating period, however it is not considered here.

The given discussion allows one to use the introduced notations together with (2.54) to obtain a matrix-vector equation relating the amplitudes of incident and diffracted modified field amplitudes in each slice of a grating layer (for details see Appendix C):

Ẽαnq = β=x,y,z N O m=-N O W αβnm Ẽβmq .
(2.58)

New matrix W replaces V in case of a corrugated grating and its components write:

W xx = ∆ -DΓ xx -DΓ xz C -1 DΓ zx , W xy = -DΓ xy -DΓ xz C -1 DΓ zy , W xz = -DΓ xz C -1 , W yx = -DΓ xy -DΓ xz C -1 DΓ zy , W yy = ∆ -DΓ yy -DΓ yz C -1 DΓ zy , W yz = -DΓ yz C -1 , W zx = -C -1 DΓ zx , W zy = -C -1 DΓ zy , W zz = I -C -1 .
(2.59)

Information about grating profile is carried by matrix Γ in form of trigonometric functions of angles ψ and ϕ. This matrix is Toeplitz relatively to the Fourier-indices m, n and writes

Γ(z) =    Γ xx Γ xy Γ xz Γ yx Γ yy Γ yz Γ zx Γ zy Γ zz    =    cos 2 ϕ sin 2 ψ sin ϕ cos ϕ sin 2 ψ cos ϕ sin ψ cos ψ sin ϕ cos ϕ sin 2 ψ sin 2 ϕ sin 2 ψ sin ϕ sin ψ cos ψ cos ϕ sin ψ cos ψ sin ϕ sin ψ cos ψ cos 2 ψ    (2.60)
Besides, additional matrices D and C that appear in (2.59) are:

D = ε ε b - ε b ε -1 (2.61) C = ε ε b -DΓ zz (2.62)
The most important difference between (2.59) and matrix V is the presence of inverse Toeplitz matrices. This means that direct replacement of V by W in Eqs. (2.47), (2.49) and (2.51) would spoil the fast numerical algorithm. To preserve the achieved advantages of the method we decompose matrix W into a product W = U(M) -1 with the following explicit form of matrices M and U (for details see Appendix C):

M =       ε b ε ε ε b 0 0 0 ε b ε ε ε b 0 0 0 ε ε b ε b ε       sin 2 ψ +    I 0 0 0 I 0 0 0 I    cos 2 ψ, (2.63) U =      ∆ x M xx + G ε ε b Γ xx G ε ε b Γ xy GΓ xz G ε ε b Γ yx ∆ y M yy + G ε ε b Γ yy GΓ yz FΓ zx FΓ zy M zz - ε b ε      , (2.64)
where

G = I - ε ε b ε b ε , F = I - ε b ε ε ε b .
(2.65)

The proposed decomposition allows one to rewrite Eqs. (2.49) and (2.51) making them free of inversions (see Appendix C):

a out = a inc + TPU(M -QRPU) -1 Qa inc (2.66) E(z q ) =     M xx 0 0 0 M yy 0 FΓ zx FΓ zy ε b ε     (M -QRPU) -1 Qa inc (2.67)
Thus, obtained equations (2.66) and (2.67) represent the required modification of (2.49) and (2.51) for the case of corrugated gratings. Eqs. (2.66) and (2.67) have more complex structure, however their formulation allows to use the proposed fast numerical algorithm with O(N log N ) time and memory resort.

Diffraction gratings in a planar structure

In the analysis of chapters 2.4 and 2.6 a diffraction grating was supposed to be placed in an isotropic homogeneous medium with the dielectric permittivity equal to the GSM basis permittivity ε b . For practical use of the method we have to generalize the obtained results for the case of arbitrary substrate and cover with permittivities ε s and ε c respectively. Let the plane interface between the basis medium and the substrate be described by amplitude reflection and transmission coefficients r for each plane harmonic with index n propagating from the inside of a grating layer. Here it will be demonstrated that multiple reflections at the mentioned interfaces can be rigorously incorporated into the method, and this does not affect the numerical complexity of the algorithm.

From the structure of matrix A one may conclude, that the changes to be introduced affect only the incident field calculation and matrices R and T. Owing the incident field amplites at a grating layer boundaries, calulation of multiple reflections is a simple exersise which leads to the following formulas: (2.68)

a e,h+ (z L ) = t (L)e,h 0 t (U )e,h 0 exp(ik z0 h) 1 -r (L)e,h 0 r (U )e,h 0 exp(2ik z0 h) a e,h+ (z U ) + r(L)e,h 0
+ (t (L)e,h 0
) 2 r (U )e,h 0 exp(2ik z0 h) 1 -r (L)e,h 0 r (U )e,h 0 exp(2ik z0 h) a e,h-(z L ), a e,h-(z U ) = r(U)e,h 0
+ (t (
Here r(L)e,h 0 and r(U)e,h 0 stand for the reflection coefficients of zero-order harmonics propagating towards a grating layer. Analogously for the waves being excited in q th sublayer we have:

a e,h+ (z q ) = t (U )e,h 0 exp [ik z0 ∆h(q -1/2)] 1 -r (L)e,h 0 r (U )e,h 0 exp(2ik z0 h) a e,h+ (z U ) + t (L)e,h 0 r (U )e,h 0 exp [ik z0 ∆h(N S + q -1/2)] 1 -r (L)e,h 0 r (U )e,h 0 exp(2ik z0 h) a e,h-(z L ), a e,h-(z q ) = t (U )e,h 0 r (L)e,h 0 exp [ik z0 ∆h(2N S -q + 1/2)] 1 -r (L)e,h 0 r (U )e,h 0 exp(2ik z0 h) a e,h+ (z U ) + t (U )e,h 0 exp [ik z0 ∆h(N S -q + 1/2)] 1 -r (L)e,h 0 r (U )e,h 0 exp(2ik z0 h) a e,h-(z L ).
(2.69) Summands 1/2 show that plane wave amplitudes and phases are measured in the middle of slices. Next, consider the propagation of harmonics from a slice with number q to a slice with number p. Presence of arbitrary substrate and cover results in the dependence of matrix R components from the direction of plane wave propagation. These components then write: and in an explicit form

R (e,h)(++) npq = ∆h θ + p-q + r (L)e,h n r (U )e,h n exp(2ik zn h) 1 -r (L)e,h n r (U )e,h n exp(2ik zn h) × exp [ik zn ∆h(p -q)] R (e,h)(+-) npq = ∆h r (L)e,h n exp [2ik zn ∆h(2N S + 1 -p -q)] 1 -r (L)e,h n r (U )e,h n exp(2ik zn h) R (e,h)(-+) npq = ∆h r (U )e,h n exp [2ik zn ∆h(p + q -1)] 1 -r (L)e,h n r (U )e,h n exp(2ik zn h) R (e,h)(--) npq = ∆h θ - p-q + r (L)e,h n r (U )e,h n exp(2ik zn h) 1 -r (L)e,h n r (U )e,h n exp(2ik zn h) × exp [-ik zn ∆h(p -q)]
(2.71)

The first and the fourth elements in (2.71) have Toeplitz structure and can be immediately multiplied by a vector via the FFT. The other two elements depend on the Fourier index sum (p + q) instead of the difference. To include them in the fast multiplication scheme one has to, first, invert a vector element numbering and change index q to index N S -q +1, then perform the FFT-based meltiplication, and finally restore the initial enumeration of the resulting vector elements. Matrix T which "gathers" diffracted harmonics and "sums" them at layer boundaries writes

T (e,h)(++) nq = exp [ik zn ∆h(N S -q + 1/2)] 1 -r (L)e,h n r (U )e,h n exp(2ik zn h) , T (e,h)(-+) nq = r (U )e,h n exp [ik zn ∆h(N S + q -1/2)] 1 -r (L)e,h n r (U )e,h n exp(2ik zn h) , T (e,h)(+-) nq = r (L)e,h n exp [ik zn ∆h(N S -q + 1/2)] 1 -r (L)e,h n r (U )e,h n exp(2ik zn h) , T (e,h)(--) nq = exp [ik zn ∆h(q -1/2)] 1 -r (L)e,h n r (U )e,h n
exp(2ik zn h) .

(2.72) Thus, Eqs. (2.68), (2.69), (2.71) and (2.72) provide the result necessary for the appli-cation of the developed method to gratings placed in an arbitrary structure. An important result of this section is the demonstration that multiple reflections at grating layer boundaries can be accounted for in a rigorous manner and can be incorporated in the fast algorithm proposed in section 2.5.

Final algorithm for the light diffraction calculation on an arbitrary profiled corrugated grating placed in a layered structure writes:

• Calculation of incident field amplitudes in each slice by Eqs. (2.68), (2.69) from given amplitudes at grating layer boundaries.

• Calculation of the Fourier-images of the dielectric and magnetic permittivities in each slice, and pre-calculation of the FFT from the obtained matrices.

• Calculation of the Fourier-images of matrix Γ (2.60) components in each slice, and pre-calculation of the FFT from the obtained matrix.

• Pre-calculation of the FFT from matrix R (2.71).

• Solution of Eq. (C.20) by the GMRES with FFT and calculation of the diffracted wave amplitudes in each slice.

• Calculation of the diffraction orders amplitudes by Eq. (2.66) and/or of the field amplitudes by (2.67).

Convergence of the numerical method

Theoretical results of the previous sections enclose the development of the fast amd memory sparing method for the light diffraction calculation on possibly very complex gratings. This section demonstrates the validity of the method by comparison of results with known refence methods. For this purpose there were chosen two methods the FMM and the Rayleigh method. The first one is rather popular and widely applicable for analysis of different diffraction structures. Additionally, both the FMM and the proposed method are Fourier methods so that one may expect them to give the same result for a given number of diffraction orders N O in the limit N S → ∞. The Rayleigh method was chosen for its perfect applicability to sinusoidal gratings [START_REF] Tishchenko | Numerical demonstration of the validity of the rayleigh hypothesis[END_REF].

To make a comparison with the FMM we considered one-and two-dimensional holographic gratings with sinusoidally changing dielectric permittivity inside the grating layer:

ε(x, y, z) =            ε c , z > z u , ε g 1 + c sin 2πx Λ x + c sin 2πy Λ y , z l ≤ z ≤ z u , ε s , z < z l , (2.73)
where c is a constant factor small enough to ensure ε ≥ 1. Grating permittivity described by Eq. (2.73) does not depend on coordinate z, and the diffraction calculation by the FMM on such grating can be performed in a particularly efficient way [START_REF] Schuster | Normal vector method for convergence improvement using the rcwa for crossed gratings[END_REF]. Fouriermatrix [ε/ε b ] mn is found analytically while calculation of [ε b /ε] mn is done numerically. For the diffraction calculation Eq. ( 2 The second benchmark was made to check the method's validity for corrugated gratings. For this purpose there was taken a sinusoidal grating with profile described by the following function: (2.75)

z s (x, y) = h 2 [sin(2πx/Λ x ) + sin(2πy/Λ y )] , (2.74)
Comparison was carried out for parameters n s = 2.5,

n c = 1, Λ x = Λ y = 1 µm, h = 0.2 µm, θ inc = ϕ inc = 30 • , λ = 0.6328 µm.
In general case the Fourier images of matrix Γ (2.60) components should be evaluated for each slice separately. However, in the particular case of sinusoidal grating there is a possibility to define them ones for all slices. Namely, the trigonometric functions of angles ψ and ϕ can be found analytically from (2.74). For 1D sinusoidal grating corresponding Fourier-integrals are found analytically. For 2D grating one can calculate the FFT of sufficiently large matrices and trace the convergence of sub-matrices of size N O1 × N O2 .

Figs. 2.6, 2.7 show the convergence for the diffraction calculation on 1D and 2D gratings respectively. As can be seen, for a fixed N O an increase of the slice number leads to a solution with some constant non-reducible error. This error demonstrates the accuracy of the infinite Fourier-sums truncation and can be reduced by an increase of N O . Examples of diffraction efficiencies calculated for the 2D sinusoidal grating are given in Tables D.3, D.4 of Appendix D. One may propose to use the power balance condition [START_REF] Petit | Electromagnetic theory of graings[END_REF] as an additional criterium of the method accuracy. However, in all simulations this condition fulfilled at least with а) б)

Λ x = Λ y = 1 µm, h = 0.2 µm, n s = 2.5, n c = 1, θ inc = ϕ inc = 30 • , λ = 0.6328 µm.
an order of magnitude better than the convergence obtained from diffraction amplitudes comparison.

Example of the calculation time dependence from the number of diffration orders is provided in Fig. 2.10 for the light diffration calculation on 2D sinusoidal grating, analogous to the previous one, however with depth 0.5 µm and N S = 250 by the developed method and the FMM. As was expected the proposed method's complexity grows approximately linearly with the number of diffraction orders. "Jumps" on the line are explained by the use of the FFT radix 2, so that the size of matrices participating in calulations doubles every time the power of 2 increases.

Conclusions

To summarize, in this chapter the generalized source method described in Chapter 1 was applied for the development of the light diffration on gratings calulcation method. The basis solution was derived for both magnetic and electric sources written in form of plane waves. From this basis solution we obtained the explicit analytical form of S-matrices of an infinitely thin grating slice. Then there was proposed an S-matrix based method for the light diffraction calculation. The method was shown to have the complexity ∼ N 3 O which is too high for complex problems solution. Therefore, the method was modified and the diffraction calculation problem was reduced to a linear system of algebraic equations with block-Toeplitz matrix. The final method has linear complexity and memory requirements with respect to calculation mesh node number. There were analized both convergence and accuracy of this method. The results presented in this chapter were published in [START_REF] Shcherbakov | Fast numerical method for modelling one-dimensional diffraction gratings[END_REF][START_REF] Shcherbakov | Fast and memory-sparing exact electromagnetic analysis of arbitrary profile 2d periodic dielectric structures[END_REF][START_REF] Shcherbakov | Fast and efficient diffraction modeling by the generalized source method[END_REF].

Chapter 3

Organic light emitting diodes with scattering layers

Light scattering calculation on nonperiodical structures

In the previous chapter we substantiated the transition from a scattering problem to a diffraction problem on gratings, and the method of exact and effective calculation of diffraction of optical radiation on 1D and 2D plane diffraction gratings was developed. Now it is necessary to return to the scattering problem in nonperiodic inhomogeneous layers. The developed approach will be applied to the solution of a problem of modeling OLEDs with a scattering layer.

Consider a nonperiodic scattering structure schematically shown in Fig. 3.1. Suppose this structure to be a homogeneous matrix containing a set of scattering particles of various size, shape, and material. Designate characteristic size of the volume filled with particles, and characteristic size of particles as Λ and d s respectively. According to the initial idea described in the introduction to Chapter 2, the developed numerical method requires transformation of geometrical parameters and coordinate functions of dielectric permittivity of a structure to the Fourier-space. Let the maximum module of a wavevector in plane XY in the considered Fourier-representation be γ max . It is necessary to choose it so that to resolve the characteristic size of particles γ max ∼ α/d s , where factor α > 1. Designating a mesh step in the reciprocal space as ∆γ, ∆γ ∼ 1/Λ we get γ max = N O ∆γ. This brings N O ∼ αΛ/d s . Choosing, for example, Λ ∼ 5µm, d s ∼ 0.5µm and α ∼ 10, one obtains the number of diffraction orders to be used N O ∼ 100.

Scattering on a nonperiodic structure described by a scattering diagram which is continuous relative to scattering angle. On the other hand the diffraction on a grating results in finite set of propagating diffraction orders. Consider a scattering volume illuminated by a plane wave propagating along axis Z. Denote a scattering amplitude as F (1) (θ) where θ is an angle between the scattering direction and axis Z (here for simplisity we consider a 2D geometry). Scattering amplitude from N identical scatteres placed along axis X equidistantly with distance Λ separating any two neighbours then writes

F (N ) (θ) = F (1) (θ) n=N/2 n=-N/2 exp(ik 0 nΛ sin θ) = 2πF (1) (θ) sin [(N + 1/2)k 0 Λ sin θ] sin k 0 Λ sin θ 2 , (3.1)
providing that we neglect the re-scattering of the radiation scattered by each volume. The last condition shows that in the limit N → ∞ the diffraction orders exactly reproduce the scattering diagramm

F (∞) (θ) = F (1) (θ) ∞ n=-∞ δ (k 0 Λ sin θ + 2πn) . (3.2)
Thus, to use solutions obtained with the developed Fourier-method to approximate corresponding solutions of scattering problems one has to decrease the influence of re-scattering.

A direct and the most simple approach to do that is to separate different scatteres to a sufficiently large distance of several wavelengths. Besides, large periods are necessary to obtain a sufficient number of diffraction orders representing a scattering diagram. Authors of [START_REF] Edee | Complex coordinate implementation in the curvilinear coordinate method: application to plane-wave diffraction by nonperiodic rough surfaces[END_REF][START_REF] Pisarenco | Aperiodic fourier modal method in contrast-feld formulation for simulation of scattering from finite structures[END_REF][START_REF] Pisarenco | Modified s-matrix algorithm for the aperiodic fourier modal method in contrast-field formulation[END_REF] additionally used perfectly matched layers (PML, [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Chew | A 3d perfectly matched medium from modified maxwell's equations with stretched coordinates[END_REF]) placed at boundaries of grating periods to simulate scattering on simple 2D dielectric bodies. Here we do not use PMLs and show that they are not necessary for scattering calculation. Now consider in detail the representation of a scattering layer within the frame of the developed method. Scattering particles with sizes comparable to the optical wavelengths in many practical applications have shape close to spherical. So further we will simulate scattering on layers with spherical particles. To define Fourier images of corresponding grating permittivities and required trigonometric functions described above, consider a simple grating containing a single spherical particle in each period. The following results can be directly generalized to the case of any number of particles.

After slicing of a grating layer period, each slice contains a circular inclusion of radius r(z) = r 2 s -z 2 (providing that the origin is placed in the center of the sphere) of a medium with ε g inside a rectangular region of permittivity ε m . Top-view of such slice is shown in Fig. 3.2. The Fourier-image of corresponding 2D function is found analytically and writes

ε mn (z) =              ε m + ∆ε πr 2 (z) Λ x Λ y , m = n = 0; ∆ε r(z)K x K y J 1 r (mK x ) 2 + (nK y ) 2 2π (mK x ) 2 + (nK y ) 2
, mn = 0.

(3.3)
where ∆ε = ε g -ε m . Calculation of the Fourier images of matrix Γ (2.60) components is Figure 3.2: Top view of a grating period slice for the grating containing a single sphere in each period. more complex. They are found for each slice separately. Placing the origin to the center of the sphere we can make the angle ψ between axis Z and normal direction to the sphere surface to be constant for each slice. Then as can be seen from (2.60) one has to consider only functions f s = sin 2 φ(x, y), f c = cos 2 φ(x, y) and f s c = sin φ(x, y) cos φ(x, y). These functions are initially defined only on the circle of raius r(z). Define them on the whole period rectangle as follows. Introduce a circular band around r(z) bounded by radii r min and r max , so that r min < r(z) < r max , and r min ≥ 0, r max < 1/2 min(Λ x , Λ y ). The period then appears to be divided into four regions, A, B, C, and D as shown in Fig. 3.2. In regions A and D take f s = f c = 0.5, f s = 0. Then introduce a "hat"-function of variable

ρ = x 2 + y 2 : f h (ρ) =                ρ r min exp - t (t-r min )(r-t) dt r r min exp - t (t-r min )(r-t) dt , r min ≤ ρ ≤ r(z), rmax ρ exp[- t (rmax-t)(t-r) ]dt rmax r exp[- t (rmax-t)(t-r) ]dt , r(z) ≤ ρ ≤ r max , (3.4)
Example of its graph is given in Fig. 3.3. Calculation of (3.4) is fast, precise and is made by Gaussian quadratures. Owing (3.4) one can write f s and f c in regions B and C as

f s (ρ) = f h (ρ) sin 2 ϕ + 1 2 [1 -f h (ρ)], f c (ρ) = f h (ρ) cos 2 ϕ + 1 2 [1 -f h (ρ)], f sc (ρ) = f h (ρ) sin ϕ cos ϕ, (3.5)
where sin ϕ = y/ρ, cos ϕ = x/ρ. Example of function f s can be seen in Fig. 3.4. Functions analogous to (3.4) are well known in the distribution theory [START_REF] Vladimirov | Generalized functions in mathematical physics[END_REF] as examples of infinitely differentiatable distributions. As one can see, definition (3.5) includes coincidense of these functions with correct values on the circle r(z), and holding of the prior trigonometric identity f s + f c = 1. Using (3.5) one obtains the necessary Fourier-matrices. The Fourier transform is made numerically by the FFT.

Scattering of a plane wave on a layer containing dielectric nanoparticles

Before application of the method one has to estimate the accuracy of the results by making a benchmark with a reference solution. One of the most well-known solutions in 3D electromagnetic wave scattering is the Mie solution [START_REF] Mie | Beiträge zur optik trüber medien, speziell kolloidaler metallösungen[END_REF] describing the scattering of a plane electromagnetic wave on a sphere. For comparison consider a grating composed of infinite number of identical dielectric spheres with refractive index 1.5 placed in a medium with refractive index 1. Take the sphere radius to be 1µm and the wavelength to be 0.6328µm. Fig. 3.5 is an illustration to the above discussion of representing of a scattering diagram by diffraction efficiencies. It can be seen that all the scattering diagram features are well reproduced. Next Fig. 3.6 demonstrates the convergence of diffraction efficiencies to the corresponding values of the differential scattering cross section (DSCS) with the increase of the diffraction order number N O for three different periods. For a fixed period the solution converges to a value with a fixed error that characterizes the amount of the re-scattered energy. This error decreases with the increase of the period and amounts to about 1% which is sufficient for comparison with most experiments.

After the benchmark consider a complex scattering layer. In an experiment the infor- mation about a scattering layer usually includes a volume or mass particle density, and particle size distribution function. The size distribution function usually can be aproximated by a Gaussian distribution with some average and dispersion. Modeling of a beam scattering on a layer as plane wave scattering implies that the diameter of the real beam is much larger than the wavelength. Besides, in the proposed method it should be much larger than the grating period. From the experimental point of view this means that the beam is scattered on a large area and an additional averaging over different groups of particles is required in the proposed method. The method calculates the diffraction fields for each wavelength and then combine them to a whole spectral response. In thick scattering layers the direct propagation of a wide spectral beam becomes incoherent on a distance exceeding several wavelengths. Thus, the typical layer thickness h is not larger than several micrometers. To treat thicker layers we divide them into N s sub-layers of thickness h, calculate the power scattering matrix components S 10 (1.14) for each sub-layer, and incrementally find the scattering matrix of the whole layer. Additionally, an averaging over different particle ensembles is included in the scattering matrix calculation with a view of getting results closer statistically to real scattering layers.

An example in Fig. 3.7 demonstrates scattering diagrams for a plane wave scattering from layers of three different thicknesses -30, 60 and 90 µm. The plane wave of wavelength 0.5 µm incidents normally to the layer. Scattering layer parameters are: ε g = 2.89, ε m = 2.67, particle diameters make a Gaussian distribution with average D s = 0.5µm and mean square δD s = 0.01µm. Grating parameters are Λ x = Λ y = 10µm and h = 3µm. One can conclude with Fig. 3.7 that the width of the zero-angle peak increases while its maximum value decreases. The scattering can be calculated via the dependence of the forward scattered power from the layer thickness h, which is given in Fig. 3.8. The last dependence is perfectly described by the exponential C exp(-κh) with attenuation coefficient κ that is found from better fitting of the curve. The dependence of κ on D s is shown in Fig. 3.9. It reveals approximately linear dependence for particle diameters that lie in range of diameters comparable to the wavelength.

One of the arguments given in Chapter 1 in favour of a rigorous method develop-

Organic light emitting diodes with scattering layers

To apply the developed method of scattering calculation for OLED simulation we have to additionally develop a general method of planar electroluminescent structures simulation. So first consider an OLED with homogeneous layers. Simulation of planar structure will be done using S-matrices described in Chapter 1. For a complete description of an elelctroluminescent planar structure Eqs. (1.15)-(1.17) should be supplemented with several relations [START_REF] Shcherbakov | Rigorous s-matrix approach to the modeling of the optical properties of oleds[END_REF]. Consider a panar structure consisting of two parts with S-matrices S U and S L separated by a dipole source layer with z-coordinate z s . Additionally, suppose that plane z = z s is not an interface separating different materials. As was written in Chapter 1 when a dipole source plane being placed in an isotropic homogeneous medium the aplitudes of the emitted field are given by (1.19). For the described structure effective amplitudes in the source plane are found from equations (according to the S-matrix definition)

a - s = S U 11 a + s + a - d , a + s = S L 22 a - s + a + d , (3.6)
which gives

a + s = 1 -S U 11 S L 22
-1 a + d + S L 22 a - d , a - s = 1 -S U 11 S L 22 -1 a - d + S U 11 a + d . (3.7)
The amplitudes of harmonics propagating outwards from the top and the bottom of the structure a - L и a + U (Fig. 3.11a), write

a + U = S U 21 a + s = S U 21 1 -S U 11 S L 22 -1 a + d + S L 22 a - d , a - L = S L 12 a - s = S L 12 1 -S U 11 S L 22 -1 a - d + S U 11 a + d . (3.8)
Next, dividing the upper part of the considered structure into two layers having matrices S U 1,2 (see Fig. 3.11b), we find amplitudes of harmonics excited in plane z = z u , z u > z s :

a + (z u ) = 1 -S U 1 22 S U 2 11 S U 2 21 a + s , a -(z u ) = 1 -S U 1 22 S U 2 11 S U 2 11 S U 1 21 a + s .
(3.9)

Analogously, in plane z = z l , z l < z s , with corresponding matrices S L1,2 (Fig. 3.11c) one finds:

a + (z l ) = 1 -S L2 11 S L1 22 S L1 22 S L2 12 a - s , a -(z l ) = 1 -S L2
11 S L1 22 S L2 12 a - s .

( Another important and necessary feature of a planar electroluminescent structure analysis is the electromagnetic power flow and loss calculation in any part of the structure. To accomplish this task we derive the Poynting vector z-projection for a given plane wave via its TE-and TM-polarization aplitudes:

S z = 1 2 ℜ{[E × H * ] z } = 1 2 ℜ(k e z ) ωµ 0 a e+ 2 -a e-2 + ℜ k h z ωε a h+ 2 -a h-2 - ℑ(k e z ) ωµ 0 ℑ a e+ a e- * -ℑ k h z ωε ℑ a e+ a e- * .
(3.11)

Note that this equation describes the energy flow for both propagating and evanscent waves in all types of layers including dielectric, lossy dielectric and metal. Analogously to Eq. (3.11) one can calculate the effective power emitted by a considered dipole layer with effective amplitudes (3.7):

P s = lim z→zs+0 S z -lim z→zs-0 S z = 1 2     ℜ(k e z ) ωµ 0 a e+ s 2 -a e+ s -a e+ d 2 + a e- s 2 -a e- s -a e- d 2 + ℜ k h z ωε a h+ s 2 -a h+ s -a h+ d 2 + a h- s 2 -a h- s -a h- d 2     - ℑ(k e z ) ωµ 0 ℑ a e+ a e- * -ℑ k h z ωε ℑ a e+ a e- * .
(3.12)

For OLED optimization it is also necessary to calculate power loss in different layers of diode structures. For a layer of thickness h bounded by planes z = z 1 и z = z 2 , z 2 -z 1 = h using (3.11) and (3.12) we get

P (z 1 , z 2 ) = P s -lim z→z 1 +0 S z -lim z→z 2 -0 S z , (3.13)
As can be seen, the equations given in this section describe only coherent wave propagation and account for all interference effects. Such description is obviously valid only when the total thickness of a structure does not exceed several wavelengths. If an OLED contains several emitting layers propagation of waves emitted by different sources should be considered incoherently (we do not study here lasing effects as, e.g., in vertically emitting cavity electroluminescent lasers [START_REF] Wilmsen | Vertical-cavity surface-emitting lasers: design, fabrication, characterization and applications[END_REF]). OLED thickness practically always allows for coherent calculation of the radiation of each separate source. However, one shoud account also for a thick substrate which impact should be included in terms of power propagation coefficients.

Presence of a substrate and corresponding subsrate-air interface results in the following consequences: all the plane waves with γ > 1 are totally reflected back to the substrate, and for γ < 1 waves are partially transmitted into the air which leads to an addition power loss. These effects are described by the reflection R a and refraction T a power coefficients for plane waves at the substrate-air interface [START_REF] Born | Principles of optics[END_REF]:

R e a = k zs -k za k zs + k za 2 R h a = k zs -ε s k za k zs + ε s k za 2 T e,h a = 1 -R e,h z , (3.14)
where k za , k z s are wavevector z-projections in the air and in the substrate, ε s is the dielectric permittivity of the substrate. Power reflection coefficient of an OLED structure is found as its S-matrix element norm R s = |S 22 | 2 . Then the total power emitted into the air P air after N reflections is related to the power emitted to the substrate P sub by

P air = P sub T a 1 -(R a R s ) N 1 -R a R s . (3.15)
The power being lost due to the re-absorbtion in OLED structure can be estimated as

P re-abs = P sub R a (1 -R s ) 1 -(R a R s ) N 1 -R a R s . (3.16)
The rest of power is trapped in the substrate:

P sub-gui = P sub (R a R s ) N . (3.17)
It should be mentioned here that power integration should be made over all wavevector projections in plane XY k x , k y . However, the integration is simplified since S-matrix components are functions of γ:

P (λ) = (2π) 2 k 4 0 ∞ -∞ ∞ -∞ P (k x , k y , λ)dk x dk y = (2π) 3 k 4 0 ∞ 0 P (γ, λ)γdγ (3.18)
Angular dependence of the emitted power in the medium with wavevector k c in the spectral region λ 1 ≤ λ ≤ λ 2 can be calulated as As a result of the derivations of this sections one may formulate a numerical method for the optical properties calculation of OLEDs with homogeneous layers:

P (θ, ϕ) = (2π) 2 k 2 c k 4
1. Setting of structure parameters: positions of layer interfaces and emitting layers, dispersion characteristics of all materials, source spectra.

2. Calculation of S-matrices for all layers and of the whole structure S-matrix by Eqs.

(1.15)-(1.18) for each wavelength λ i , i = 0, . . . , N λ -1, and plane wavecector projection γ i , i = 0, . . . , N k -1, from 0 to some γ max with step ∆γ = γ max /N k . Benchmark for the last method was established by comparison with experimental results obtained for a test OLED [START_REF] Shcherbakov | Rigorous s-matrix approach to the modeling of the optical properties of oleds[END_REF]. Structure of this OLED is shown in Fig. 3.12 and includes a green emitting layer, electron and hole blocking and transporting layers, metal cathode and transparent ITO anode. Figs. 3.13 and 3.14 show the comparison between numerical and experimental results for the spectral emitted power and OLED color coordinates angular dependence. As can be seen the correspondence between experimental and simulation data is quite good. Additionally, Fig. Next consider the same green OLED and introduce in it a scattering layer between the anode and the substrate as shown in Fig. 3.12. Let the scattering layer to consist of a homogeneous matrix with embedded spherical particles with refractive index constrast equal to 0.1. Take the particle diameter to be 0.5µm and particle volume density 0.1. Fig. 3.16 shows the simulated spectral dependence of the OLED external efficiency for several thicknesses h sc of the scattering layer. It is seen that introduction of a thin scattering layer allows to increase the efficiency from fractions of a percent to about 6% depending on the wavelength. With increase of the wavelenght the efficiency also increases. This feature can be explained on the basis of the modal analysis of the OLEd structure. Fig. 3.17 shows spatial field distributions for zero TE and TM modes. TM mode is plasmonic one and is concentrated at the metal-dielectric interface between cathode and electron transporting layer. Due to its high localization and low propagation length it can hardly be scattered. On the contrary, the field distribution of the TE mode is much broader and shifts towards the scattering layer location with the increase of the wavelenght. Besides, the increase of the wavelength leads to decrease of the TE mode propagation constant, thus, exponentially increasing the efficiency of the TE mode scattering (see Fig. 3.10). Fig. 3.16 shows that scattering layers of thicknesses up to dozens of microns do scatter trapped radiation and increase the OLED efficiency by several percent of the absolute Figure 3.13: Spectral emitted power for the device shown in Fig. 3.12 for measurement and simulation with and without a half-ball lens. value. Increase of the scattering layer thickness leads to an increase of scattering and corresponding decrease of the efficiency which is demonstrated in Fig. 3.19.

Provided simulation examples allow to make the following conclusions. First, mamximum gain that can be obtained by introduction of a scattering layer is achieved for rather thin layers from several microns to several dozens of microns. Second, to achieve the external efficiency increase by 5-10% one has to carry out an optimization of both OLED and scattering layer parameters. Such optimization makes sense only for specific device configurations and materials and was not done in this work.

Conclusions

This chapter describes a method of a nonperiodic scattering problem solution via simulation of the light diffraction on corresponding periodized structure. The proposed numerical benchmark consists in comparison of the diffraction efficiencies calculated for 2D grating of dielectric spheres with the DSCS given by the Mie solution, and shows the possibility to get about 1% accuracy for the range of problems under consideration. Next we developed a numerical S-matrix based method for the planar electroluminescent structures simulation. Validity of the method was demonstrated by comparison with experimental results. Finally both methods for scattering calculation and OLED simulation were joined, and there was demonstrated an increase of the OLED efficiency due to scattering layer inclusion with accompanying discussion of mechanism of this effect.

Conclusion

In this final chapter we list and discuss the main results of the thesis. Chapter 1 is devoted to the review of the works related to problems solved in the thesis. First, the range of methods capable to rigorously solve diffraction and scattering problems in plane layered spatially nonhomogeneous structures is listed, and shortages of these methods are specified. The representation used in Fourier-methods being close to the approach, developed in this work, is described in more details. Special attention is paid to the S-matrix technique used in Chapter 3 for modeling of organic light-emitting diodes. Secondly, section 1.3 describes the generalized source method previously proposed by the supervisor of the thesis author in [START_REF] Tishchenko | A generalized source method for wave propagation[END_REF][START_REF] Tishchenko | Generalized source method: new possibilities for waveguide and grating problems[END_REF]. The GSM became a cornerstone of the theoretical model developed in Chapter 2. The last section of the review contains the description of results on OLED optical properties simulation. There we describe a problem of the light outcoupling from planar OLED structures and provide arguments in favor of use of scattering layers. On the basis of the review it is shown that a rigorous modeling of optical properties of OLED with scattering layers is quite difficult with existing methods, and there is a need in new fast and memory sparing method.

Chapter 2 is devoted to the application of the GSM to the light diffraction on gratings problem. In the first paragraph there are given considerations substantiating the relation between scattering and diffraction on gratings problems relying on the requirement of fast matrix-vector multiplications in the Fourier space. Then, on the basis of the GSM we analytically derive S-matrix components for infinitely thin slices of diffraction gratings. These formulas allowed to formulate the method of the light diffraction calculation. However, as was shown the method has too high numerical complexity to be concurrent and useful O(N 3 O ), with N O being the number of points in the Fourier space. To reduce the numerical complexity the method was reformulated so as to calculate wave amplitudes in all slices of a grating at once by solving a system of linear equations. This system is derived for both holographic and corrugated gratings. There was proposed a numerical algorithm relying on the GMRES where matrix-vector multiplications are found by the FFT. This was shown to be possible due to the Toepltiz structure of the linear system matrix. The resulting complexity of the method was estimated to be O(N O N S ), where N S is the slice number. This result represents a prominent step in the numerical methods development for the problem of the light diffraction on gratings. In the end of the Chapter 2 we demonstrate the analysis of the method convergence and accuracy.

In the first section of Chapter 3 we discuss the back transition from the solution of a diffraction on a grating problem to the corresponding soution of the scattering problem. Validity of the proposed approach is demonstrated by the benchmark with the Mie solution. We showed that range of problems under consideration can be solved with about 1% accuracy or better. Next we developed a numerical S-matrix based method for optical simulation of the planar electroluminescent structures. Validity of the method was demonstrated by comparison with experimental results. Finally both methods for scattering calculation and OLED simulation were joined, and there was demonstrated an increase of the OLED efficiency due to scattering layer inclusion with accompanying discussion of mechanism of this effect.

To conclude we point out one more time that the main result of the thesis is the munerical method for the diffraction on gratings calculation with linear numerical complexity relative to the number of mesh nodes. Calculation of OLED optical properties is one of possible aplications of the method. Obviously it perfectly suites for simulation of diffractive structures with fully determined geometry, diffractive optical elements. In addition, the derived equations allow to consider not only electric sources but also magnetic, or, in other word, variations in the magnetic permittivity. This gives means to directly include perfectly matched layers into models and, probably, improve solution of non-periodic problems.

a e± = E ± (k × ẑ) |k × ẑ| = E ± x k y γ -E ± y k x γ (A.
3)

a h± = H ± (k × ẑ) |k × ẑ| = H ± x k y γ -H ± y k x γ (A.4)
where γ = k 2 x + k 2 y . Using the Maxwell's equations one gets

a e± = ± ωµ b k x 2γk z H x ± ωµ b k y 2γk z H y - ωµ b 2γ H z (A.5) a h± = ∓ ωε b k x 2γk z E x ∓ ωε b k y 2γk z E y + ωε b 2γ E z (A.6)
On the basis of (A. Additionally, introduce matrix Q, composed of (A.7), (A.8): Here the first sign "±" is specified by incident harmonics, and the second by diffracted harmonics. Thus, the derived Eqs. (B.2)-(B.10) represent general analytical relations for S-matrix components of an initely thin slice of a corrugated grating.

E H = Q E Q H a =

Appendix C Derivation of formulas describing the light diffraction on corrugated gratings

This appendix describes in detail the derivation of equations of section 2.3, Chapter 2, describing the electromagnetic wave diffraction on corrugated gratings. The introduced local coordinate system allows to explicitly write Eq. (2.54): Note that application of the Fourier transform to Eq. (C.6) is correct since all the above dericaltions preserve the coorect order of factors in multiplications as well as difference between the normal and tangential fields and sources. As noticed in section 2.6 the derived matrix W contains inverse Toeplitz sub-matrices C -1 and cannot be applied in the proposed fast numerical algorithm. To solve this problem consider the structure of matrix C -1 :

   j n j ψ j φ    = -iωε b       ε b ε -1 -I 0 0 0 ε ε b -I 0 0 0 ε ε b -I          E n E ψ E φ    (C.
C -1 = ε ε b sin 2 ψ + ε b ε -1 cos 2 ψ -1 = M -1 xx ε b ε = M -1 yy ε b ε = ε b ε M -1 zz , (C.7)
where matrix M is defined by (2.63). Eqs. (C.7) and (2.65) lead to

C -1 D = -FM -1 xx = -FM -1 yy , DC -1 = -GM -1
xx .

(C.8)

These relations allow to rewrite components of W as follows together with (2.63) and (2.64). Now the derived decomposition (C.18) should be used to reformulate the linear system of equations.

W xx = ∆ x -DΓ xx -DΓ xz C -1 DΓ zx = ∆ x M xx + G ε ε b Γ xx M -1 xx , (C.
To prove (2.66) we substitute (C.18) into (2.47) and (2.49). The first one defining amplitudes of diffracted harmonics in the grating layer is rewritten as [START_REF] Wasey | Effects of dipole orientation and birefringence on the optical emission from thin films[END_REF] Notice that this equation includes non-rectangular matrices. We do not specify the size of unit matrices as they are uniquely defined by other ones. In correspondence with the general GSM formulation presented in section 1.3, unknown wave amplitudes inside the grating layer are found from self-consistent equation:

Q E a inc = Q E (I -RPWQ E )a = (I -Q E RPUM -1 )Q E a = (M -Q E RPU)M -1 Q E a. (C.
a = a inc + RPU(M -Q E RPU) -1 Q E a inc (C.20)
A direct substitution justifies that (C.20) meets (C.19). Finally, Eq. (C.20) allows to calculate the diffracted harmonics amplitudes at a grating layer boundaries:

a out (z U , z L ) = a inc (z U , z L ) + TPWQ E a = a inc (z U , z L ) + TPUM -1 Q E a inc + TPUM -1 Q E RPU(M -Q E RPU) -1 Q E a inc = a inc (z U , z L ) + TPU I + M -1 Q E RPU(I -M -1 Q E RPU) -1 M -1 Q E a inc = a inc (z U , z L ) + TPU(I -M -1 Q E RPU) -1 M -1 Q E a inc = a inc (z U , z L ) + TPU(M -Q E RPU) -1 Q E a inc
E q =    I 0 0 0 I 0 C -1 DΓ zx C -1 DΓ zy C -1    Q E a =    I 0 0 0 I 0 C -1 DΓ zx C -1 DΓ zy C -1    M(M -Q E RPU) -1 Q E a inc =     M xx 0 0 0 M yy 0 -FΓ zx -FΓ zy ε b ε     (M -Q E RPU) -1 Q E a inc

Figure 1 . 1 :

 11 Figure 1.1: Plane layered structure with Cartesian coordinates.

  methods. Now define S-and T-matrices by means of the introduced notations. Consider plane waves propagaing in the substrate and in the cover. Let their amplitudes be a ± s and a ± c respectively (Fig. 1.1). Then, a T-matrix relates wave amplitudes in the substrate and in the cover and writes a +

2 × 2 ,Figure 1 . 2 :L 11 - 1 S L 10 SL 11 - 1 S U 01 SL 11 - 1 S U 10 S 11

 21211110111011111011 Figure 1.2: To definitions of a) a plane interface S-matrix, b) a homogeneous layer Smatrix and c) an S-matrix multiplication rule.

Figure 1 . 3 :

 13 Figure 1.3: Example of a 2D periodic sinusoidal intefrace separating two media within a plane layer.

Figure 1 . 4 :

 14 Figure 1.4: Examples of OLEDs: a) conventional OLED and b) an OLED with a scattering layer.

Figure 1 . 5 :

 15 Figure 1.5: Methods of the OLED efficiency increase: a) conventional OLED, b) OLED with aerogel, c) corrugated OLED, d) nanostructured substrate, e) OLED with microlenses, and f) OLED with mesastructure [96].

. 25 )

 25 According to the GSM, generalized currents are proportional to products of fields by permittivity modifications(1.34),(1.35). First, consider an index grating described by continuous functions ε(r), µ(r) of coordinates x, y. Relation between the Fourier-components of the fields and generalized currents follows from (1.34), (1.35), (2.22), (2.23), and writesj x,ym = -iωε b ε x,y ε b mn -δ mn Ẽx,yn , j zm = -iωε b δ mn -ε b ε z mn Ẽzn ,(2.26)f x,ym = iωµ b µ x,y µ b mn -δ mn Hx,yn , f zm = iωµ b δ mn -µ b µ z mnHzn .(2.27)Now one can write out explicit formulas for S-matrix components of an infinitely thin slice of a plane grating

2 .

 2 Calculation of matrices containing the Fourier harmonics of dielectric and magnetic permittivities ε(x, y, z p ), µ(x, y, z p ) in each slice p = 0 . . . N L -1 with z p = z (L) + (p + 1/2)∆h;3. Calculation of matrices containing theFourier harmonics of thrigonometric functions of angles defining normal directions at curved interfaces separating different media inside a grating layer (see Appendnices B, C); 4. Calculation of S-matrices of each slice using Eqs. (2.30)-(2.33) or (B.2)-(B.10); 5. Calculation of the whole grating S-matrix by means of Eq. (1.18).

  and, generally, this operation is made by O(N 3 O ) multiplications. Since Eqs.(1.18) are used each time when a new slice is added, the net numerical complexity of the method appears to be O(N3

Figure 2 . 1 :Figure 2 . 2 :

 2122 Figure 2.1: Convergence of the S-matrix-based diffraction calculation method for a rectangular grating. a) Convergence as a function of N S for different N O . b) Use of polynomial interpolation to improve the convergence.

Figure 2 . 3 :

 23 Figure 2.3: Local Cartesian coordinates n, ψ, ϕ at an interface between different media, with axes directions defined by angles ψ и ϕ.

  z0 h) a e,h-(z L ).

  .49) is used. Examples of results obtained by the developed method and the FMM are given in Tables D.1, D.2 of Appendix D for the following parameters: n g = n s = 2.5, n c = 1, c = 1, Λ x = Λ y = 1 µm, h = 0.5 µm, θ inc = ϕ inc = 30 • , λ = 0.6328 µm. Figs. 2.4, 2.5 demonstarte the convergence of the method with the increase of the slice number N S and comparison with the result obtained by the FMM. Number of the Fourier harmonics N O is not specified since within the limits of N O from 10 to 150 for 1D grating and from N O = N OX N OY = 5 × 5 = 25 to N O = N OX N OY = 50 × 50 = 1600 for 2D grating the obtained dependecies were quite similar.

Figure 2 . 4 :

 24 Figure 2.4: Convergence (relative error) and comparison with the FMM (absolute error) versus the inverce slice number for the light diffraction calculation on 1D holographic sinusoidal grating. Parameters of the problem are: Λ x = 1 µm, h = 0.5 µm, c = 0.1, n s = n g = 2.5, n c = 1, θ inc = 30 • , λ = 0.6328 µm.

Figure 2 . 5 :

 25 Figure 2.5: Convergence (relative error) and comparison with the FMM (absolute error) versus the inverce slice number for the light diffraction calculation on 2D holographic sinusoidal grating. Parameters of the problem are: Λ x = Λ y = 1 µm, h = 0.5 µm, c = 0.1, n s = n g = 2.5, n c = 1, θ inc = ϕ inc = 30 • , λ = 0.6328 µm.

Figure 2 . 6 :

 26 Figure 2.6: Convergence with the slice number increase for the light diffraction calculation on 1D corrugated sinusoidal grating. Parameters of the problem are: Λ x = 1 µm, h = 0.2 µm, n s = 2.5, n c = 1, θ inc = 30 • , λ = 0.6328 µm.

Figure 2 . 7 :

 27 Figure 2.7: Convergence with the slice number increase for the light diffraction calculation on 2D corrugated sinusoidal grating. Parameters of the problem are: Λ x = Λ y = 1 µm, h = 0.2 µm, n s = 2.5, n c = 1, θ inc = ϕ inc = 30 • , λ = 0.6328 µm.

Figure 2 . 8 :

 28 Figure 2.8: Maximum absolute difference between calculated diffraction amplitudes by the GSM and the Rayleigh method from the slice number for the light diffraction calculation on 1D corrugated sinusoidal grating for a) TE-and b) TM-polarization. Parameters of the problem are: Λ x = 1 µm, h = 0.2 µm, n s = 2.5, n c = 1, θ inc = 30 • , λ = 0.6328 µm.

Figure 2 . 9 :

 29 Figure 2.9: Maximum absolute difference between calculated diffraction amplitudes by the GSM and the Rayleigh method from the slice number for the light diffraction calculation on 2D corrugated sinusoidal grating. Parameters of the problem are: Λ x = Λ y = 1 µm, h = 0.2 µm, n s = 2.5, n c = 1, θ inc = ϕ inc = 30 • , λ = 0.6328 µm.

Figure 2 . 10 :

 210 Figure 2.10: Calculation time dependence from the number of diffration orders for the light diffration calculation on 2D sinusoidal grating for the GSM and the FMM.

Figure 3 . 1 :

 31 Figure 3.1: Calculation problems of scattering on a single object and diffraction on the same periodized object.

Figure 3 . 3 :Figure 3 . 4 :

 3334 Figure 3.3: Example of a graph of function f h , defined by (3.4), for r min = 1 and r max = 2.

Figure 3 . 5 :

 35 Figure 3.5: Normalized differential scattering cross section provided by the Mie solution of a plane wave scattering on a sphere, and normalized diffraction efficiencies for a plane wave diffraction on a grating composed of the same spheres for a) TE polarization and b) TM polarization of the incident wave.

Figure 3 . 6 :

 36 Figure 3.6: Convergence to the Mie solution.

Figure 3 . 7 :

 37 Figure 3.7: Scattering diagrams for a plane wave scattering from layers of three different thicknesses -30, 60 and 90 µm.

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Dependence of the forward scattered power from the layer thickness h.

Figure 3 . 10 :

 310 Figure 3.10: Dependence of the evanescent wave scattering cross section from the in-plane wave projection of this wave. Solid line correspond to function σ(γ) = C 1 exp(-C 2 γ) with coeeficients C 1 = 12.5, C 2 = 13.4.

. 10 )Figure 3 . 11 :

 10311 Figure 3.11: Calculation of plane wave amplitudes inside a planar structure: a) effective source amplitudes b) amplitudes in the upper part of a structure relative to the source layer c) amplitudes in the lower part of a structure relative to the source layer.

3 .

 3 Calculation of wave amplitudes at all interfaces by Eqs. (3.7)-(3.10) for all wavelengths λ i and in-plane wavevector projections γ i . 4. Calculation of power loss in each layer, effective power for each source, and power emitted to the subsutrate and to the air by Eqs. (3.11)-(3.13).

  3.15 demonstrates the power loss analysis carried out by Eqs. (3.11)-(3.13).

Figure 3 .

 3 Figure 3.12: OLED example used for comparison of experimental and simulation data.

Figure 3 . 14 :Figure 3 . 15 :

 314315 Figure 3.14: Measured and simulated color coordinates for the device shown in Fig. 3.12.

Figure 3 . 16 :

 316 Figure 3.16: Spectral external efficiency for the OLED, shown in Fig. 3.12, with and without a scattering layer.

Figure 3 . 17 :

 317 Figure 3.17: Modal field spatial distribution of zero TE and TM modes for the OLED, shown in Fig. 3.12. Vertical lines indicate interfaces between different OLED layers, and 0 corresponds to the cathode-air interface.

Figure 3 . 18 :

 318 Figure 3.18: Spectral dependence of the zero TE mode propagation constant for the OLED, shown in Fig. 3.12.

Figure 3 . 19 :

 319 Figure 3.19: Dependence of the external efficiency of the OLED, shown in Fig. 3.12, with introduced scattering layer from the thickness of this layer. Horizontal line shows the efficiency of the device without scattering layer.

  3)-(A.6) it is possible to derive the inverse transform which writes

1 ) 1 5 )

 115 Using transformations (2.55),(2.56) rewrite this equation in the initial coordinates XY Z iωε b (∆ x I αβ -DΓ αβ ) matrices ∆ x , D and Γ αβ are defined by (2.29), (2.61) and (2.60) respectively. Toderive the relation between source j and modified field Ẽ write out separately zcomponent of (C.2):j z = -iωε b [-DΓ zx E x -DΓ zy E y + (∆ x -DΓ zz )E z ] . (C.3)Then use the definition of the modified field (2.22) and (2.62) to getẼz = E z -j z iωε b = -DΓ zx E x -DΓ zy E y + ε ε b -DΓ zz E z = -DΓ zx Ẽx -DΓ zy Ẽy + CE z . (C.4)This enables one to explicitly specify relation between the real and the modified fields insode the grating layer:DΓ zx C -1 DΓ zy C -1Then substituting (C.5) into (C.2) we get equation providing the explicit form of matrixW E (2.59): x -DΓ xx -DΓ xz C -1 DΓ zx -DΓ xy -DΓ xz C -1 DΓ zy -DΓ xz C -1 -DΓ xy -DΓ xz C -1 DΓ zy ∆ y -DΓ yy -DΓ yz C -1 DΓ zy -DΓ yz C -1 -C -1 DΓ zx -C -1 DΓ zy I -C -1

9 )W

 9 xy = -DΓ xy -DΓ xz C -1 DΓ zy = G ε ε b Γ xy M -1 yy , (C.10) W xz = -DC -1 Γ xz = GΓ xz M -1 zz , (C.11) W yx = -DΓ yx -DΓ yz C -1 DΓ zx = G ε ε b Γ yx M -1 xx , (C.12) W yy = ∆ y -DΓ yy -DΓ yz C -1 DΓ zy = ∆ y M yy + G ε ε b Γ yy M -1 yy , (C.13) W yz = -DC -1 Γ yz = GΓ yz M -1 zz , (C.14) W zx = -C -1 DΓ zx = FΓ zx M -1 xx , (C.15)W zy = -C -1 DΓ zy = FΓ zy M -1 yy , (C.[START_REF] Lukosz | Theory of optical-environment-dependent spontaneous emission rates for emitters in thin layers[END_REF])W zz = I -C -1 = M zz -

(C. 21 )

 21 Using (C.5), (C.20), (C.7) and (C.8) one finds amplitudes of the electric field in each slice:

(C. 22 )Table D. 2 :R - 1 ,T 3 :R - 1 ,T - 1 4 :

 22213114 Diffraction efficiencies calculated by the GSM for the diffraction of a plane wave on a holographic sinusoidal grating with the following parameters of the problem:n s = n g = 2.5, n c = 1, c = 0.1, Λ x = Λ y = 1 µm, h = 0.5 µm, θ inc = ϕ inc = 30 • . Wavelength λ = 0.6328 µm. T E → T E T E → T M T M → T E T M → T M 0.00002923 0.00001101 0.00005734 T -1,-1 0.00434517 0.00179887 0.00203993 0.00471443 T -1,00.00130657 0.05517556 0.06046262 0.00091777 Diffraction efficiencies calculated by the Rayleigh method for the diffraction of a plane wave on a corrugated sinusoidal grating with the following parameters of the problem:n s = 2.5, n c = 1, Λ x = Λ y = 1 µm, h = 0.2 µm, θ inc = ϕ inc = 30 • . Wavelength λ = 0.6328 µm. T E → T E T E → T M T M → T E T M → T M 0.00012129 0.00012129 0.06583007 R 0,1 0.00769935 0.00939271 0.00344653 0.00265069 T -1,-1 0.00585463 0.00352645 0.00412859 0.01043321 Diffraction efficiencies calculated by the GSM for the diffraction of a plane wave on a corrugated sinusoidal grating with the following parameters of the problem n s = 2.5, n c = 1, Λ x = Λ y = 1 µm, h = 0.2 µm, θ inc = ϕ inc = 30 • . Wavelength λ = 0.6328 µm. T E → T E T E → T M T M → T E T M → T M R -1,-1 0.00570901 0.00165777 0.00170973 0.00538269 R -

  2 . This enables Eqs. (2.18) and (2.19) to be rewritten as (here we refer to the introduced modified field (2.22), (2.23))

  zm γ m k xm ∆ mn k yn -k xm ∆ mn k xn -k mx Ω xxmn k ny + k mx Ω xymn k nx -k my Ω yxmn k ny + k my Ω yymn k ny 1 γ n -ωε b γ m Σ zxmn k yn -Σ zymn k xn zm γ m k xm ∆ mn k xn + k ym ∆ mn k yn -k mx Ω xxmn k nx -k mx Ω xymn k ny -k my Ω yxmn k nx + k my Ω yymn k ny k zn γ n + 1 k zm γ m k xm Υ xzmn + k ym Υ yzmn γ n zm γ m k xm ∆ mn k xn + k ym ∆ mn k yn -k mx Ω xxmn k nx -k mx Ω xymn k ny -k my Ω yxmn k nx + k my Ω yymn k ny k zn γ n + 1 k zm γ m k xm Υ xzmn + k ym Υ yzmn γ n zm γ m k xm ∆ mn k xn + k ym ∆ mn k yn -k mx Ω xxmn k nx -k mx Ω xymn k ny -k my Ω yxmn k nx + k my Ω yymn k ny k zn γ n -1 k zm γ m k xm Υ xzmn + k ym Υ yzmn γ n zm γ m k xm ∆ mn k xn + k ym ∆ mn k yn -k mx Ω xxmn k nx -k mx Ω xymn k ny -k my Ω yxmn k nx + k my Ω yymn k ny k zn γ n -1 k zm γ m k xm Υ xzmn + k ym Υ yzmn γ n DΓ αβ + DΓ αz C -1 Γ zβ D, Υ αβ = DΓ αβ C -1 , Σ αβ = C -1 Γ αβ D (B.11) in accordingly to (2.61), (2.62) and (2.60). Also, consider separately an important case of collinear TM diffraction on one-dimensional grating. Taking ϕ = 0 and k xm=0 in (2.55), (2.56) one obtains DΓ yy -DΓ yz C -1 DΓ zy mn (±k zn ) + 1 ±k zm DΓ yz C -1 mn k yn + 1 k ym C -1 DΓ zy mn k nz + 1 k ym I -C -1 mn k yn .

	S h-e+ mn	= S h-e-mn	=	ωε b k 1 γ m	,
											(B.6)
	S h+h+ mn		= k + 1 γ m	1	Σ zxmn k xn + Σ zymn k yn	k zn γ m	+	1 γ m	I -C -1	mn γ n ,	(B.7)
	S h+h-mn		= -1 γ m k +		1 Σ zxmn k xn + Σ zymn k yn	k zn γ m	+	1 γ m	Qa I -C -1	mn γ n ,	(A.9) (B.8)
	S h+h-mn		= -1 γ m k +		1 Σ zxmn k xn + Σ zymn k yn	k zn γ m	+	1 γ m	I -C -1	mn γ n ,	(B.9)
	S h+h-mn	= k + 1 γ m	1	Σ zxmn k xn + Σ zymn k yn	k zn γ m	+	1 γ m	I -C -1	mn γ n ,	(B.10)
	where
	Ω αβ = S h±h± mn = 1 k zm ∆ -

(B.12)

Appendix A

Plane wave polarization

In this application we introduce notations defining amplitudes of TE-and TM-polarized waves via the electric and the magnetic field compponent amplitudes, and corresponding invers tranformations. Solutions of the Maxwell's equations in an isotropic homogeneous medium with permittivities ε b and µ b can be written in form of plane waves:

Here the sign "±" distinguishes waves propagating in the positive and negative directions relative to axis Z, and the net field amplitude is a sum of amplitudes with signs "+" and "-". Wavevector projections obey the condition [START_REF] Born | Principles of optics[END_REF]

In general case when k z can be a complex value, waves with 0 ≤ arg k z < π are denoted by "+", and waves with π ≤ arg k z < 2π by sign "-". Define amplitudes of TE a e± and TM a h± waves as Appendix B

S-matrices of corrugated gratings

This appendix describes the derivation of components of a profiled grating S-matrix. For this purpose use Eq. ( 2

The last equaition defines the required S-matrix (2.28) which components explicitly write

Appendix D

Tables of diffraction efficiencies

Table D.1: Diffraction efficiencies calculated by the FMM for the diffraction of a plane wave on a holographic sinusoidal grating with the following parameters of the problem: