D. Gottlieb, S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications (Society of Industrial and Applied Mathematics, Philadelphia, Pa., 1989).

H. P. Herzig, M. T. Gale, H. W. Lehmann, R. Morf, “Diffractive components: computer generated elements,” in Perspectives for Parallel Optical Interconnects, P. H. Lalanne, P. Chavel, eds. (Springer-Verlag, Berlin, 1993), Chap. 5, pp. 71–107.

[CrossRef]

S. E. Sandström, G. Tayeb, R. Petit, “Some comments on the use of exact eigenfunctions for lossy lamellar gratings,” Internal rep. (Laboratoire d’Optique et Electromagnétique, Centre de Saint-Jérôme, 13397 Marseille Cedex 20, France, 1992);cf. also S. E. Sandström, G. Tayeb, R. Petit, “Lossy multistep lamellar gratings in conical diffraction mountings: an exact eigenfunction solution,” J. Electromagn. Waves Appl. 7, 631–649 (1993).

[CrossRef]

Its convergence properties appear to be unsatisfactory for gratings with very large coupling amplitudes or for cases in which spectrally very narrow resonances occur. To my knowledge, no attempt has been made thus far to use it for metallic gratings.

R. Petit, ed., Electromagnetic Theory of Gratings, Vol. 22 of Topics in Current Physics (Springer-Verlag, Berlin, 1980).

[CrossRef]

The research cited in Ref. 11 is based on Fourier approximations for the refractive index (cf. Ref. 2), which lead to unsatisfactory convergence behavior. This result is avoided in the work of Li,10 who uses exact eigenfunctions for lamellar gratings. However, Li does not attempt to identify the source of instability or to explain why his R-matrix propagation procedure works. Furthermore, instabilities are still encountered “whose real cause is still unknown.”10

C. Heine, R. Morf, “Submicrometer gratings for solar energy applications,” Appl. Opt.34 (to be published).

[PubMed]

Actually, it is easy to see that expanding a discontinuous function in terms of analytic functions must give rise to poor convergence: everyone is familiar with the Fourier transformation of the square-wave function, whose Fourier coefficients an decrease as 1/n. Similarly, we can see that a function that is k− 1 times continuously differentiable but whose kth derivative is discontinuous would have Fourier amplitudes that decay as 1/nk+1. Consequently, if a function is infinitely many times continuously differentiable, its Fourier amplitudes will decay faster than any finite power of the order n.

The matrix is, in general, complex and nonsymmetric. However, in contrast to the case of Fourier approximation,2 for the case of H polarization no generalized eigenvalue problem results. I use software routine eigcc(from the IMSL Numerical Libraries distributed by Visual Numerics Inc., Houston, TX) for computing eigenvalues and eigenvectors, which takes order M3 operations, where M is the rank of the matrix, i.e., the number of polynomial coefficients specifying the eigenfunctions.

M. Abramowitz, I. A. Stegun, Handbook of Mathematical Functions, Vol. 55 of NBS Applied Mathematics Series (National Bureau of Standards, U.S. Government Printing Office, Washington, D.C., 1972), pp. 437–438.

Care must be taken, though, that for low-order approximations the order Ml does not drop below 2, in which case the differential equation is no longer solved and only boundary conditions are transmitted across the layer. In our calculations we use Ml ≥ 4.

The usual method, employed by Botten et al.3 consists of solving the adjoint problem (replacing the dielectric constant by its complex conjugate). Our method, based on polynomial expansion, involves the calculation of both the left-hand and the right-hand eigenvectors of the algebraic eigenvalue problem.

L. Collatz, The Numerical Treatment of Differential Equations (Springer-Verlag, Berlin, 1960), pp. 49–50.

The slow decay of the matrix elements is, of course, the reason that so many eigenvectors must be included in the calculation. Although they may not transmit information from one interface to the next because of their fast decay in space (for high orders), these elements still affect the distribution of amplitudes at the interface at which they are excited. Thus only a stable code can accurately take their effect into account.

It may be useful to consider a simple example that illustrates this question very clearly. Assume that we wish to compute two different sums: S(N)=∑n=1N(1/n) and T(N)=∑n=1N[1/(n+1)], whose difference is 1 − 1/(N + 1). A direct calculation of S(10,000) − T(10,000) (in single-precision arithmetic) leads to an error of 3 × 10−6if we carry out the summation starting at n = 1 and increase n, whereas the reverse order leads to an error more than 25 times smaller. By contrast, if we carry out the summation by successively adding 1/n and subtracting 1/(n + 1), there is an exact cancellation of error. Thus, although this calculation is clearly a stable procedure, its precision depends crucially on the order of steps.

There is no a priori reason to keep the number of eigenmodes the same in each grating layer. Indeed, we would expect to have to use more eigenfunctions in grating layers with high refractive indices exhibiting a greater number of propagating modes and less rapid spatial decay for the evanescent waves. Equation (60) allows us to have the freedom to define the number of modes individually. However, this freedom is not utilized in the present study and is the subject of future research.

J. H. Wilkinson, C. Reinsch, Linear Algebra (Springer-Verlag, Berlin, 1971), pp. 70–92.In these algorithms column pivoting, which is compatible with the band shape, is used to improve numerical stability.

J. H. Wilkinson, The Algebraic Eigenvalue Problem (Clarendon, Oxford, 1965), pp. 233–236.The triangularization can also be accomplished by standard matrix factorization (LU factorization, with L and U denoting lower and upper triangular matrices, respectively), based on the Crout algorithm or on Gaussian elimination.

A trivial type of instability occurs if one of the eigenvalues χn becomes exactly zero. Then the equations relating the field at interfaces between grating layers are singular. This result occurs for special isolated wavelengths. We can avoid this problem by using wavelengths in the vicinity of these isolated wavelengths for which none of the χn becomes too small.

If the eigenvalues and the eigenvectors are already known, the computation of matrices that do not depend on grating depth requires a total of 9N3 + O(N2) complex multiplications and additions per grating layer for E polarization, whereas for H polarization 4N3/3 additional operations are required for the computation of matrix N(j)[Eqs. (47) and (A78)]. These numbers include the triangularization of matrices C3(j) by Householder transformations and the associated modifications of matrices C4(j) and −I, which are necessary to reduce the number of subdiagonal elements from 2N to N.

In most of the numerical research relating to diffraction theory no detailed analysis of the operation count of numerical procedures is made. For the comparison of different algorithms operation count is a universally useful criterion, more so than CPU time, which depends on the particular implementation and on the special features of the processor. Lacking this information for other algorithms, we can use CPU times as a comparison: Li specifies the CPU time on a SPARCstation 2 for a calculation with N = 60 eigenmodes for a 10-layer system as approximately 200 s, or 20 s/layer (for E or H polarization).10 In the present method, with the same machine type, the computation times for each grating layer are as follows (the numbers in parentheses are for a Digital Equipment Corporation Alpha-AXP 3000-500): 24 (3.6) s for the algebraic eigenvalue problem of dimension 3N = 180, 8 (1.2) s for the calculation of Nmost important eigenvectors, 6 (0.7) s for the computation of the grating depth hj-independent matrices, whereas the (hj-dependent) solution of Eq. (60) amounts to 2.5 (0.32) s/layer, all by double-precision arithmetic. When single-precision arithmetic is used the time for calculation of eigenvalues and eigenvectors on the SPARC is reduced from a total of 32 to approximately 20 s, whereas the remaining steps are reduced by only roughly 10%. On the Alpha these computation times are cut down by approximately a factor of 2 if single-precision arithmetic is used. These times are valid for FORTRAN programs that use standard compilers without using blaslibraries or other optimized codes.

I. S. Gradsteyn, I. M. Ryzhik, Table of Integrals, Series, and Products (Academic, New York, 1965), p. 836.

The justification for the approximation used in Appendix B is as follows. Assume that we introduce a very thin layer of vacuum between successive grating layers, in which the eigenmodes would be only the Fourier waves. Then our insertion of a complete set of plane waves becomes physically justified. If this layer is of the order of a few monolayers of typical dielectrics, i.e., in the nanometer range, its effect becomes essentially invisible; however, it will introduce some cutoff in the number of Fourier waves that contribute significantly.