From our own experience, we are able to warn the nonspecialist, that a rough numerical application of the integral formalism cannot furnish reliable numerical results.  Prof. Daniel Maystre, "Rigorous vector theories of diffraction gratings", in Progress in Optics XXI, E. Wolf, ed., Elsevier Science Publishers B.V., 1984, pp. 167.


An accurate numerical solution of the problem becomes more complicated when we encounter at least one of the following conditions: several interfaces between different transmitting and reflecting materials; distance between the interfaces that is comparable or smaller than typical dimensions of the structure; these interfaces are not equal in their length; there are inclusions inside the layers; the layers are anisotropic; the structure is nonperiodic but long; the form of the incident wave front is not plane; the incident beam is oblique (conical diffraction); the flux in the finite light beam has a Gaussian distribution; the source of light is distributed; the grating is not plane; the grating periodicity is two or threedimensional; and so on. Accurate methods have become particularly favored in recent years. In theory, almost all of them enable one to solve the problems of diffraction by complicated periodic structures, even though they usually represent a plane wave approximation over a limited range of parameters. ^{ }The efficiency of diffraction gratings used at present in spectroscopy and photonics essentially depends on the borders profiles and their depths, on layered materials and their thickness, on a range of wavelengthtoperiod ratio, and on a choice of working wavelengths and/or orders. Concurrent optimization of all these parameters cannot be realized in practice without numerical modeling. We review here only methods which have been developed since the mid1960's, and we intend to find a solution without making any guesses as to mathematical and physical formulation of the problem or to its solution algorithm. Unlike the approximate theories, restrictions are imposed only on the stages of numerical implementation and computation. They all are based on Maxwell's equations, exact boundary conditions, radiation conditions and are well known as exact (or rigorous) electromagnetic theories. ^{}^{}^{}^{}^{42}^{}
The most widely used accurate methods of grating diffraction analysis are as follows: the integral method; ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{18}^{}^{,}^{}^{}^{}^{}^{42}^{}^{}^{}^{}^{}^{}^{44}^{} the classical differential method; ^{}^{}^{}^{}^{42}^{}^{,}^{}^{}^{}^{}^{45}^{} the modal method; ^{}^{}^{}^{}^{46}^{}^{}^{}^{}^{}^{}^{48}^{} the coupledwave approach, which combines an exact coupledwave analysis and Fourier modal method; ^{}^{}^{}^{}^{45}^{}^{,}^{}^{}^{}^{}^{49}^{}^{}^{}^{}^{}^{}^{51}^{} the coordinate transformation approach, which includes the socalled "C" method and conformal mapping method; ^{}^{}^{}^{}^{45}^{}^{,}^{}^{}^{}^{}^{52}^{}^{}^{}^{}^{}^{}^{54}^{} the finiteelement or finitedifference approaches, which include the finiteelement method, the finitedifference time domain method, and the boundary element method; ^{}^{}^{}^{}^{55}^{}^{}^{}^{}^{}^{}^{58}^{} and the multipole method. ^{}^{}^{}^{}^{59}^{}^{,}^{}^{}^{}^{}^{60}^{}
Within the large variety of accurate theoretical approaches and their modifications which are used to calculate the diffraction gratings efficiency, the integral equation method seems to hold the lead. Nowadays this accurate numerical method ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{18}^{}^{,}^{}^{}^{}^{}^{42}^{}^{}^{}^{}^{}^{}^{44}^{}^{ }is universally recognized as one of the most developed and flexible techniques for solving diffraction problems. ^{}^{}^{}^{}^{1}^{}^{,}^{}^{}^{}^{}^{32}^{}^{,}^{}^{}^{}^{}^{45}^{}^{,}^{}^{}^{}^{}^{61}^{} Historically, this method enabled one for the first time to accurately solve vector problems of light scattering by optical diffraction gratings and to achieve an excellent agreement with experimental data. ^{}^{}^{}^{}^{42}^{} This was due to the high accuracy and good convergence of the method, ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{44}^{}^{,}^{}^{}^{}^{}^{62}^{} especially for TM polarization planes. In fact many wellknown optical problems of diffraction by periodic and nonperiodic structures were originally handled with the theory of integral equations. Thus far this method has been one of the most powerful tools for analysis of diffraction by almost any type of grating over the entire wavelength range ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{62}^{}^{,}^{}^{}^{}^{}^{63}^{} despite very intensive development of other exact methods. ^{}^{}^{}^{}^{32}^{}^{,}^{}^{}^{}^{}^{45}^{}^{,}^{}^{}^{}^{}^{53}^{}^{,}^{}^{}^{}^{}^{59}^{} In many important cases, the integral method is the only acceptable approach from the practical standpoint for performing research and making accurate predictions of the efficiency characteristics ^{}^{}^{}^{}^{1}^{}^{,}^{}^{}^{}^{}^{15}^{}^{,}^{}^{}^{}^{}^{61}^{}^{,}^{}^{}^{}^{}^{63}^{}^{, }^{}^{}^{}^{}^{64}^{}. At the same time, presentday progress of the integral method, as well as of other exact approaches, goes hand in hand with development of more diversified numerical algorithms.
The integral method is an approach which allows us in a rigorous manner to reduce a problem of diffraction by grating to solving a linear boundary integral equation or a system of such equations.^{}^{}^{}^{}^{42}^{} In general the integral approach, as well as the similar finiteelement method, implies twodimensional integration. However, in actual practice, a onedimensional curvilinear integration easily reduced to ordinary integrals is used. Then the linear integral equations so obtained are reduced to a system of linear algebraic equations by the collocation method ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{44}^{} or by Galerkin's method ^{}^{}^{}^{}^{65}^{}. In our realization, we use a rather simple but robust and universal technique, the modification of socalled classical Nyström's method ^{}^{}^{}^{}^{66}^{}. The process of numerical solution of integral equations is based on discretization with piecewise constant weighting functions. The principal parameter, in which the convergence is estimated, is the number N of discretization (collocation) points on each boundary. The collocation points and the quadrature nodes can be chosen independently of one another or selected at the same locations. The latter choice (our program's default option in most cases) requires a standard regularization of integrals. ^{}^{}^{}^{}^{44}^{} It is worth noting that the regularization is used even at corner nodes of a nonsmooth boundary. ^{}^{}^{}^{}^{43}^{} For relatively shallow profiles, the nodes can be uniformly arranged along the xcoordinate (grating periodicity). But a uniform distribution with respect to the arc length, as shown in, ^{}^{}^{}^{}^{44}^{}^{,}^{}^{}^{}^{}^{62}^{} is more universal and makes it possible to treat, for example, lamellar, severely asymmetrical, or even nonfunction profiles by the integral method without any additional effort at the user's end. In a narrow sense, the modified integral method (MIM) is a discretization method which also specifies a summation rule for the kernel functions. In the simplest case, the corresponding series is truncated symmetrically at the lower summation index  P and the upper index + P, where P is an integer defined by P ~ gN. The “truncation ratio” g is optimized at small values of N and is kept constant as N increases. It has been found ^{}^{}^{}^{}^{62}^{} that g = 1 /2 is a reasonably good option for many practical calculations. Quadratures in our codes are performed by the rectangular rule with the following singleterm corrections: for the Green function  we take into account its logarithmic singularity; for normal derivatives of the Green function  we account for the profile curvature. ^{}^{}^{}^{}^{42}^{} For gratings with smooth boundaries, this method yields the overall error estimate O( N^{3}) for diffraction amplitudes and efficiencies in both polarization planes. However, the above simple truncation rule for the Green series is inadequate to match such accuracy of the discretization. An efficacious remedy is provided by Aitken's acceleration procedure ^{}^{}^{}^{}^{67}^{} applied to the truncated Green series.
To date the integral method ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{44}^{}^{ }has been perhaps the only accurate method which enables one to rather easily study the efficiencies of gratings with real groove profiles in any spectral range. It is due to the nature of the method itself when the path of integration for the surface current density is coincident with real surfaces of the grating layer boundaries. This is particularly true for the MIM, ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{62}^{} in which groove profiles are represented not in a distorted form in terms of the Fourier expansions (about smoothing edge effects see ^{}^{}^{}^{}^{42}^{}), but in a more correct form by means of discretization points which coincide with points of the real groove profile. Hence we can scrupulously take into account all jumps of the border profile functions and their first and second derivatives. Moreover, this strategy makes it possible to calculate the efficiencies of gratings with real multilayer border profiles which have fine structures, including random micro and nanoroughness. In the recent past, a generalization of the integral method has also been proposed for arbitrary rough gratings and incident beams (Gaussian beam, in particular) ^{}^{}^{}^{}^{5}^{}^{,}^{}^{}^{}^{}^{8}^{}. ^{ }The most general integral theory enables ^{}^{}^{}^{}^{32}^{} us to deal with almost all kinds of gratings problems, including such difficult cases as: any kind of echelles ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{14}^{}^{,}^{}^{}^{}^{}^{15}^{}^{,}^{}^{}^{}^{}^{64}^{}^{,}^{}^{}^{}^{}^{61}^{}^{,}^{}^{}^{}^{}^{68}^{} (used in orders up to 1431, or even higher ^{}^{}^{}^{}^{15}^{}); bulk and multilayer gratings with real groove profiles in the Xray  EUV range; ^{}^{}^{}^{}^{1}^{}^{,}^{}^{}^{}^{}^{2}^{}^{,}^{}^{}^{}^{}^{4}^{}^{,}^{}^{}^{}^{}^{5}^{}^{,}^{}^{}^{}^{}^{7}^{}^{,}^{}^{}^{}^{}^{13}^{}^{,}^{}^{}^{}^{}^{29}^{}^{}^{}^{}^{}^{}^{31}^{} very deep gratings with arbitrary profiles ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{44}^{}, especially at high conductivity and in the TM polarization; ^{}^{}^{}^{}^{62}^{} etc. Among the disadvantages of this method are some mathematical complicacy and various "particularities" of its numerical realization. Specifically, an approximation of the Green functions and their normal derivatives is applied, with resulting accuracy sufficient for common practice and with satisfactory computation time. ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{18}^{}^{,}^{}^{}^{}^{}^{44}^{}^{,}^{}^{}^{}^{}^{62}^{} In addition, the integral method may at times be difficult to adapt to heterogeneous or anisotropic media.
Differential Method
Modal Method
The classical modal (also referred to as the characteristic wave or characteristicmode approach) method (MM) ^{}^{}^{}^{}^{46}^{}^{}^{}^{}^{}^{}^{48}^{} offers an accurate formulation of the problem without any supplementary assumptions. In this method, Maxwell's equations are solved in a closed form in terms of field expansions for each rectangular (steplike or lamellar) layer in closed form, with the solutions in different layers connected by corresponding boundary conditions. Compared to the coupledwave method, the modal method refers to merely alternative mathematical representations of the field inside the grating. In their expanded rigorous forms both formulations are completely equivalent. However, the modal approach is more natural, since the total field inside the grating is expressed as a weighted sum taken over all possible characteristic waves ("modes"). Each individual mode satisfies Maxwell's equations with corresponding boundary conditions and propagates through the periodic medium unchanged. It thus becomes possible to evaluate the electromagnetic field inside the grooves to a greater accuracy. The modal method allows one to easily evaluate highconducting lamellar profiles. ^{}^{}^{}^{}^{32}^{}^{,}^{}^{}^{}^{}^{47}^{} The method can be used not only for rectangulargroove gratings. Various profiles can be analyzed by partitioning into thin layers parallel to the surface. ^{}^{}^{}^{}^{45}^{}^{,}^{}^{}^{}^{}^{47}^{}^{ }This approach has several restrictions ^{}^{}^{}^{}^{32}^{}^{,}^{}^{}^{}^{}^{62}^{} and is most effective in the case of lamellar gratings.
CoupledWave Method
The rigorous coupledwave analysis (RCWA) ^{}^{}^{}^{}^{49}^{}^{}^{}^{}^{}^{}^{51}^{} (sometimes called the method of Moharam and Gaylord or the Fourier modal method  FMM) is contiguous with the classical differential and modal methods. In each separate rectangular slice, the solution of the ordinary differential equations with constant coefficients can be found as a sum over "modes". The mode coefficients are determined by matching the field expansions at the slice boundaries. This is a direct application of the classical differential formalism to lamellar profiles. A solitary coupled wave does not satisfy Maxwell's equations. As in the differential method, this approach requires Fourier representations of the field components and the permittivity at both sides of the corrugated region and hence cannot deal with highly reflecting surfaces. Despite the provision of significantly improved convergence using the fast Fourier factorization (FFF) ^{}^{}^{}^{}^{45}^{} and similar techniques, ^{}^{}^{}^{}^{50}^{}^{,}^{}^{}^{}^{}^{51}^{} both methods still bring about weak convergence for high conducting nearinfrared gratings in the TM polarization. ^{}^{}^{}^{}^{69}^{} As in the modal method, nonlamellar profiles can be represented in the form of several rectangular steps and treated using the RCWA. However, this brings about numerical instability ^{}^{}^{}^{}^{32}^{} and, in particular, a loss of precision, which produces a need to use some special techniques, like S and Rmatrix algorithms, ^{}^{}^{}^{}^{70}^{} orthonormalization, and others, ^{}^{}^{}^{}^{49}^{} as is the case for the classical differential formalism. The method is most suitable for dielectric gratings, for which a small number of Fourier harmonics are required, and a stepwise approximation does not give rise to large errors in the efficiencies.
Coordinate Transformation Method
The idea behind the coordinate transformation method (also known as Chandezon method) or closely similar conformal mapping method resides in introducing a new coordinate system that not only maps the corrugated surfaces to planar surfaces, which simplifies the boundary conditions, but also transforms Maxwell's equations in the Fourier space into a matrix eigenvalue problem. ^{}^{}^{}^{}^{53}^{} Thus, Maxwell's equations in the new curvilinear nonorthogonal coordinate system are recast into equations with variable coefficients which admit of a straightforward numerical solution of the particular grating problem. This method can deal with deep gratings regardless of polarization and the refractive index, as well as with multilayer nonconformal ^{}^{}^{}^{}^{52}^{} and inhomogeneous anisotropic gratings. ^{}^{}^{}^{}^{54}^{} The coordinate transformation method does not cope effectively with some limiting cases of small wavelengthtoperiod ratios and especially of high orders (echelles). ^{}^{}^{}^{}^{32}^{} In addition, weak convergence is observed for real border profiles which have abrupt jumps of the profile function derivatives. ^{}^{}^{}^{}^{32}^{}^{, }^{}^{}^{}^{}^{62}^{}
Finite Element Group of Methods
Rigorous methods of analysis, which are the most universal but, at the same time, the most cumbersome for developing corresponding algorithms and getting numerical results, are being used for some kinds of diffraction gratings or nonperiodic complex structures when there is no other alternative. The wellknown versions of the appropriate theories are widely used for solving Maxwell's equations by finiteelement (FE) or finitedifference (FD) methods. ^{}^{}^{}^{}^{56}^{} The finiteelement method, the method of finitedifference time domain (FDTD), and the boundary element method (BEM) must be included in the list of these approaches. All three are conceptually different realizations of the same theory and thus they have much in common. In these methods the entire field is represented as a sum of elementary functions over the mesh cells number. The explicit choice of such mesh functions allows one to analytically integrate corresponding differential equations, reducing the problem to a system of linear algebraic equations for the unknown amplitudes. The major distinction from the aforementioned rigorous approaches resides in the fact that the twodimensional Maxwell PDEs are solved numerically on the sampled grid. However, such direct twodimensional integration is known for instability and demands for a great deal of computational resources. Not long ago these methods have been adjusted for solving the problem of the light dispersion by gratings. ^{}^{}^{}^{}^{55}^{}^{}^{}^{}^{}^{}^{58}^{} Being rather complicated and unadapted to solution of specific diffraction grating problems, these approaches are also characterized by numerical inaccuracy which is due to the twodimensional meshes used. ^{}^{}^{}^{}^{55}^{} Moreover, bulk linear systems must be used to solve some typical problems with the required accuracy. ^{}^{}^{}^{}^{57}^{}
Multipole Method
The multipole approach (sometimes called the method of fictitious sources or the multiple multipole (MMP) method) is a rigorous theory which allows one to eliminate the singularities in the integral equations by neglecting them when they are positioned above or below the profile (usually above and below concurrently). ^{}^{}^{}^{}^{59}^{}^{,}^{}^{}^{}^{}^{60}^{} Since the region in which the diffracted field is evaluated is separated from the fiction source region, singularities that come from both the source and the viewing point approaching each other do not exist. By using this approach, the diffraction problem can be reduced to evaluating (by least square fit) both the amplitudes of the electromagnetic field radiated by the fictitious sources (poles of the scattering matrix) and those of the incident field. The method of fictitious sources seems at first glance to avoid the disadvantages of both the differential method (passing "across" the profile) and the integral method (the complexity of isolating singularities). However, a price is paid by the user of the code: there is no reliable general algorithm for choosing the position of the fictitious surfaces as well as the position and the type of the fictitious sources. ^{}^{}^{}^{}^{32}^{} Despite some peculiarities in utilization, the multipole method is a suitable tool for solving scattering problems on complex structures such as photonic crystals and waveguide modes. ^{}^{}^{}^{}^{59}^{}^{,}^{}^{}^{}^{}^{60}^{}
The interested reader can find more about the existing electromagnetic theories and their implementation in ^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{32}^{}^{,}^{}^{}^{}^{}^{42}^{}^{,}^{}^{}^{}^{}^{45}^{}^{,}^{}^{}^{}^{}^{53}^{} or in recent publications devoted to the matter at this website or on the Internet.
Nowadays several dozen analytical, semianalytical and numerical methods exist to predict grating efficiency behavior. But there are only several commercially available tools for an exact modeling of diffraction by gratings. Some of them are presented here.
Table. Commercial codes for grating efficiency modeling
Name

Method

Web Address

Founder

Year

Price

PCGrate™

MIM


Dr. L. I. Goray

1989

€1000€20,000

MaX1™

MMP & FD


Prof. C. Hafner

1990

£650

GSolver™

RCWA & MM


Dr. D. Fluckiger

1994

$495$1,995.

DiPoG™

GFEM


Dr. J. Elschner, Dr. R. Hinder, Dr. A. Rathsfeld, Dr. G. Schmidt

2003

€1,000€20,000

GratingMOD™

RCWA


Dr. R. Scarmozzino

2004

~$10,000

GDCalc™

RCWA


Mr. K. C. Johnson

2005

$550

They are all based both on the aforementioned rigorous approaches and on their different numerical implementation. However, before buying them, it is well to know how the code covers parameters of your gratings and how the method converges in your immediate practical applications. Under these circumstances, the program will meet the requirements you imposed on the design and the tapping of such a sophisticated device as any grating is. This can result in great savings in time and money both in the manufacturing environment and at the user's laboratory.^{}^{}^{}^{}^{71}^{} The result obtained by using the code should improve the quality of gratings and other instruments, and also hasten the cycle of their development. True enough, the reliability of the results thus obtained is somewhat maintained by the totality of experience with all these programs. A detailed comparison of various commercial programs and scientific codes in terms of efficiency measurements is presented in ^{}^{}^{}^{}^{5}^{}^{,}^{}^{}^{}^{}^{6}^{}^{,}^{}^{}^{}^{}^{9}^{}^{,}^{}^{}^{}^{}^{15}^{}^{,}^{}^{}^{}^{}^{44}^{}^{, }^{}^{}^{}^{}^{48}^{}^{,}^{}^{}^{}^{}^{53}^{}^{,}^{}^{}^{}^{}^{61}^{}^{}^{}^{}^{}^{}^{63}^{}.

