From our own experience, we are able to warn the non-specialist, 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. 1-67.

Any optical system can be correctly described by the corresponding Maxwell equations with appropriate boundary conditions which take into consideration both characteristic properties of the radiation source and discontinuities of the electromagnetic field at the interfaces between media components in the system of interest. Then one should remember that physical characteristics (and, in particular, the refractive index) are continuous inside each medium, while at the interfaces sharp changes of these characteristics may occur. Approximations and methods used for simplified solutions of Maxwell’s equations are examined comprehensively in41. If characteristic dimensions of the diffraction element significantly exceed the wavelength of the incident light, then, as a rule (but not always!), these types of structures are easily analyzed by approximate numerical methods such as the scalar Kirchhoff-type simulation. Dimensions of modern diffraction structures are comparable or less than the wavelengths of propagating waves. In this case, the corresponding boundary value problem cannot be solved by any of the approximate or simplified approaches. Even if an approximate asymptotic method yields good results for specific parameters of the problem (in comparison with the accurate method or with an experiment), it can give totally erroneous predictions when these parameters are changed slightly. In other words, it is impossible to evaluate an error introduced by any approximate method for the vector diffraction problem prior to the accurate solution of the problem being obtained. A clear understanding of this facilitates rapid development of the accurate numerical procedures for solving problems of diffraction by periodic structures.42

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 non-periodic 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 three-dimensional; 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 wavelength-to-period 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 mid-1960’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 coupled-wave approach, which combines an exact coupled-wave analysis and Fourier modal method;45,49–51 the coordinate transformation approach, which includes the so-called “C” method and conformal mapping method;45,52–54 the finite-element or finite-difference approaches, which include the finite-element method, the finite-difference time domain method, and the boundary element method;55–58 and the multi-pole 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 method9,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 well-known optical problems of diffraction by periodic and non-periodic 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 range9,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 characteristics1,15,61,63, 64. At the same time, present-day 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 finite-element method, implies two-dimensional integration. However, in actual practice, a one-dimensional 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 method9,44 or by Galerkin’s method65. In our realization, we use a rather simple but robust and universal technique, the so-called classical Nyström collocation method66. The process of numerical solution of integral equations is based on collocation with piecewise constant basis functions. The principal parameter, in which the convergence is estimated, is the number N of 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 non-smooth boundary.43 For relatively shallow profiles, the nodes can be uniformly arranged along the x-coordinate (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 non-function 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 collocation method which also specifies a summation rule for the Green functions and their normal derivatives. 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 ~ kN. The “truncation ratio” k is optimized at small values of N and is kept constant as N increases. It has been found62 that k = 1/2 is a reasonably good option for many practical calculations. Quadratures in our codes are performed by the rectangular rule with the following single-term 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 collocation. An efficacious remedy is provided by Aitken’s acceleration procedure67 applied to the truncated Green series.

To date the integral method9,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 see42), but in a more correct form by means of collocation 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 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-roughness. 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 enables32 us to deal with almost all kinds of gratings problems, including such difficult cases as: any kind of echelles9,14,15,64,61,68 (used in orders up to 1431, or even higher15); bulk and multilayer gratings with real groove profiles in the X-ray – EUV range;1,2,4,5,7,13,29–31 very deep gratings with arbitrary profiles9,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.

BRIEF DESCRIPTION OF OTHER

The classical differential method is an accurate method by which we can reduce the problem of diffraction by grating to resolution of partial differential equations with appropriate boundary conditions.42 The projection method is used to transform these two-dimensional equations into a set of coupled ordinary differential equations. The truncated equations are integrated by well known numerical algorithms.45 The natural basis that is usually used is the Fourier representations of both the field and the permittivity inside the corrugated region. The main feature of the differential method is that it “works” across the profile using respective Fourier transforms for constant horizontal lines. Due to specific boundary conditions in the TM polarization, the vertical derivative of the field has a jump when crossing the corrugated surface. This requires a large number of Fourier components in the expansion of the field for highly conducting surfaces despite recent significant advances in convergence acceleration in the TM plane.69 The differential method is amenable for treating many (but not all) gratings efficiency problems.32,45,69

The classical modal (also referred to as the characteristic wave or characteristic-mode 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 (step-like or lamellar) layer in closed form, with the solutions in different layers connected by corresponding boundary conditions. Compared to the coupled-wave 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 high-conducting lamellar profiles.32,47 The method can be used not only for rectangular-groove gratings. Various profiles can be analyzed by partitioning into thin layers parallel to the surface.45,47 This approach has several restrictions32,62 and is most effective in the case of lamellar gratings.

The rigorous coupled-wave 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 near-infrared gratings in the TM polarization.69 As in the modal method, non-lamellar profiles can be represented in the form of several rectangular steps and treated using the RCWA. However, this brings about numerical instability32 and, in particular, a loss of precision, which produces a need to use some special techniques, like S- and R-matrix 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.

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 non-orthogonal 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 non-conformal52 and inhomogeneous anisotropic gratings.54 The coordinate transformation method does not cope effectively with some limiting cases of small wavelength-to-period 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

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 non-periodic complex structures when there is no other alternative. The well-known versions of the appropriate theories are widely used for solving Maxwell’s equations by finite-element (FE) or finite-difference (FD) methods.56 The finite-element method, the method of finite-difference 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 two-dimensional Maxwell PDEs are solved numerically on the sampled grid. However, such direct two-dimensional 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 two-dimensional meshes used.55 Moreover, bulk linear systems must be used to solve some typical problems with the required accuracy.57

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 in9,32,42,45,53 or in recent publications devoted to the matter at this website or on the Internet.

COMMERCIAL CODES

Nowadays several dozen or even hundred numerical computer packages exist to solve Maxwell’s equations and predict grating efficiency behavior. But there are only several commercially available specialized tools for an exact modelling of diffraction by spectral gratings. Some of them are presented here.

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 in5,6,9,15,44, 48,53,61–63.

Name Method Web Address Founder Year Price
Comsol Multiphysics™ FEM/BEM comsol.com Dr. S. Littmarck & Mr. F. Saeidi 1986 $/€4,000–$/€40,000
PCGrate™ MIM/BIM pcgrate.com
iigrate.com		 Prof. L. I. Goray 1989 $/€2,000–$/€35,000
GSolver™ RCWA & MM gsolver.com Dr. D. Fluckiger 1994 $795
DiPoG® BEM wias-berlin.de Dr. J. Elschner et al. 2003 ~ $10,000
GratingMOD™ RCWA rsoftdesign.com Dr. R. Scarmozzino 2004 ~ $10,000