ACCURATE
ELECTROMAGNETIC THEORIES
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 in^{41}. 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 Kirchhofftype 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 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,4244} the classical differential method;^{42,45} the modal method;^{4648} the coupledwave approach, which combines an exact coupledwave
analysis and Fourier modal method;^{45,4951} the coordinate transformation approach, which includes the
socalled "C" method and conformal mapping method;^{45,5254} the finiteelement or finitedifference approaches, which include
the finiteelement method, the finitedifference time domain method, and the boundary
element method;^{5558} 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,4244
}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.
INTEGRAL METHOD
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 socalled classical Nyström collocation method^{66}. 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 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
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 found^{62} 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 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 collocation. 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 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 microroughness. 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,2931} 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.
BRIEF DESCRIPTION OF OTHER METHODS
Differential Method
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 twodimensional
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}
Modal Method
The classical modal (also referred to as the characteristic wave or
characteristicmode approach) method (MM)^{4648} 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)^{4951} (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.^{5558} 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.
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.
Table. Commercial codes for grating efficiency
modeling
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,6163}.
