Abstract

A new numerical mathematical method is presented for stack of lamellar gratings. It is based on expansion of the electromagnetic field in terms of the eigenfunctions of the Helmholtz equation. These eigenfunctions are in turn expanded in terms of sets of polynomial basis functions. It is shown that, for arbitrary polarization and for both dielectrics and lossy metals, the eigenvalues and the eigenvectors and, consequently, the spectral location of resonances converge at an exponential rate with increasing dimension of the polynomial basis. For the solution of the boundary-value problem physical arguments are used to derive a new algorithm that is of high numerical accuracy and is inherently stable. Single-precision arithmetic is sufficient, even for the calculation of strong resonances.

© 1995 Optical Society of America

Full Article  |  PDF Article

References

  • View by:
  • |
  • |
  • |

  1. R. Petit, ed., Electromagnetic Theory of Gratings, Vol. 22 of Topics in Current Physics (Springer-Verlag, Berlin, 1980).
    [CrossRef]
  2. K. Knop, “Rigorous diffraction theory for transmission gratings with deep rectangular grooves,” J. Opt. Soc. Am. 68, 1206–1210 (1978).
    [CrossRef]
  3. L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
    [CrossRef]
  4. L. C. Botten, M. S. Craig, R. C. McPhedran, “Highly conducting lamellar gratings,” Opt. Acta 28, 1103–1106 (1981).
    [CrossRef]
  5. D. Gottlieb, S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications (Society of Industrial and Applied Mathematics, Philadelphia, Pa., 1989).
  6. 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]
  7. 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]
  8. D. M. Pai, K. A. Awada, “Analysis of dielectric gratings of arbitrary profiles and thicknesses,” J. Opt. Soc. Am. A 8, 755–762 (1992);cf. also M. Neviére, E. Popov, “Analysis of dielectric gratings of arbitrary profiles and thicknesses: comment,” J. Opt. Soc. Am. A 9, 2095–2096 (1992).
    [CrossRef]
  9. 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.
  10. L. Li, “Multilayer modal method for diffraction gratings of arbitrary depth and permittivity,” J. Opt. Soc. Am. A 10, 2581–2591 (1993).
    [CrossRef]
  11. N. Château, J. P. Hugonin, “Algorithm for the rigorous coupled-wave analysis of grating diffraction,” J. Opt. Soc. Am. A 11, 1321–1331 (1994).
    [CrossRef]
  12. 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
  13. C. Heine, R. Morf, “Submicrometer gratings for solar energy applications,” Appl. Opt.34 (to be published).
    [PubMed]
  14. 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.
  15. 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.
  16. 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.
  17. 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.
  18. 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.
  19. L. Collatz, The Numerical Treatment of Differential Equations (Springer-Verlag, Berlin, 1960), pp. 49–50.
  20. 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.
  21. 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.
  22. 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.
  23. 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.
  24. 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.
  25. 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.
  26. 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.
  27. 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.
  28. I. S. Gradsteyn, I. M. Ryzhik, Table of Integrals, Series, and Products (Academic, New York, 1965), p. 836.
  29. 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.

1994 (1)

1993 (1)

1992 (1)

1981 (2)

L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
[CrossRef]

L. C. Botten, M. S. Craig, R. C. McPhedran, “Highly conducting lamellar gratings,” Opt. Acta 28, 1103–1106 (1981).
[CrossRef]

1978 (1)

Abramowitz, M.

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.

Adams, J. L.

L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
[CrossRef]

Andrewartha, J. R.

L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
[CrossRef]

Awada, K. A.

Botten, L. C.

L. C. Botten, M. S. Craig, R. C. McPhedran, “Highly conducting lamellar gratings,” Opt. Acta 28, 1103–1106 (1981).
[CrossRef]

L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
[CrossRef]

Château, N.

Collatz, L.

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

Craig, M. S.

L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
[CrossRef]

L. C. Botten, M. S. Craig, R. C. McPhedran, “Highly conducting lamellar gratings,” Opt. Acta 28, 1103–1106 (1981).
[CrossRef]

Gale, M. T.

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]

Gottlieb, D.

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

Gradsteyn, I. S.

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

Heine, C.

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

Herzig, H. P.

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]

Hugonin, J. P.

Knop, K.

Lehmann, H. W.

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]

Li, L.

McPhedran, R. C.

L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
[CrossRef]

L. C. Botten, M. S. Craig, R. C. McPhedran, “Highly conducting lamellar gratings,” Opt. Acta 28, 1103–1106 (1981).
[CrossRef]

Morf, R.

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]

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

Orszag, S. A.

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

Pai, D. M.

Petit, R.

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]

Reinsch, C.

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.

Ryzhik, I. M.

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

Sandström, S. E.

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]

Stegun, I. A.

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.

Tayeb, G.

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]

Wilkinson, J. H.

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.

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. Opt. Soc. Am. (1)

J. Opt. Soc. Am. A (3)

Opt. Acta (2)

L. C. Botten, M. S. Craig, R. C. McPhedran, J. L. Adams, J. R. Andrewartha, “The dielectric lamellar grating,” Opt. Acta 28, 413–428 (1981);“The finitely conducting lamellar grating,” 28, 1087–1102 (1981).
[CrossRef]

L. C. Botten, M. S. Craig, R. C. McPhedran, “Highly conducting lamellar gratings,” Opt. Acta 28, 1103–1106 (1981).
[CrossRef]

Other (23)

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.

Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (18)

Fig. 1
Fig. 1

Geometry of the diffraction problem.

Fig. 2
Fig. 2

Illustration of propagating modes am′+(i+1) and am′(i+1) generated at the top interface. These modes couple information directly across the grating. By contrast, the decaying modes am+(i) and am(i) are characterized by their exponential decay away from the interface of generation. They do not couple directly across the grating structure, but they do modify the boundary condition at the interface at which they are generated. Note that the higher the spatial resolution of the representation of the electromagnetic field, the more and the faster decaying modes are generated.

Fig. 3
Fig. 3

Convergence of the third eigenvalue for the geometry (cf. Fig. 1) with c = 0.4, d = 1, with dielectric constants 2 = out= 1, 1= 25 for E polarization with normal incidence. Note the slow decay of the error for Fourier approximation, as 1/M2, where M specifies the order. By contrast, both Chebyshev and Legendre polynomial expansions lead to exponential convergence of the error. Note also that for very low order the Fourier result is superior.

Fig. 4
Fig. 4

Convergence of the third eigenvalue for H polarization. Here the error decays as 1/M for Fourier approximation. By contrast, as for E polarization (cf. Fig. 3), both Chebyshev and Legendre polynomial expansions lead to exponential convergence of the error. Note that, in contrast to the case of E polarization, the error of the Fourier expansion is larger even for very small order M.

Fig. 5
Fig. 5

Total transmission for the geometrical and the optical parameters of Fig. 3 at normal incidence for a grating height h = 0.1: E polarization in Legendre representation.

Fig. 6
Fig. 6

Total transmission for the geometrical and the optical parameters of Fig. 3 at normal incidence for a grating height h = 0.1: H polarization in Legendre representation.

Fig. 7
Fig. 7

Convergence of transmission for E polarization at the edge of resonance λ = 1.25d (cf Fig. 5).

Fig. 8
Fig. 8

Convergence of transmission for H polarization at the edge of resonance λ = 1.23d (cf. Fig. 6).

Fig. 9
Fig. 9

Resonance at λ = 1.14d of total transmission in E polarization. Fourier expansion (cf. Fig. 5).

Fig. 10
Fig. 10

Resonance at λ = 1.13d of total transmission in H polarization. Fourier expansion (cf. Fig. 6).

Fig. 11
Fig. 11

Resonance at λ = 1.14 d of total transmission in E polarization. Legendre expansion (cf. Fig. 9).

Fig. 12
Fig. 12

Resonance at λ = 1.13d of total transmission in H polarization. Legendre expansion (cf. Fig. 10).

Fig. 13
Fig. 13

Deep rectangular grating (Fig. 7 of Ref. 8). Transmission in E polarization of orders 0 and −1. Legendre expansion with 11 and 23 modes. Note the nonzero transmission of order 0 near Θ = 30°, in contrast with the results obtained in Ref. 8.

Fig. 14
Fig. 14

Same as Fig. 13, with an expanded scale for angle of incidence. Notice the minimum of transmission order −1 at 18°, which is significantly lower than the result obtained in Ref. 8.

Fig. 15
Fig. 15

Zero-order transmission for metallic grating. The geometric parameters are c = 0.1, d = 1, h = 0.1. The dielectric constants are all unity, with the exception of 1 = (1.8 + 7.12i)2. The angle of incidence is Θ =5°. Note the Wood anomaly at λ = 0.913. Legendre expansion.

Fig. 16
Fig. 16

Same as Fig. 15, with an expanded wavelength scale. Note the excellent convergence of Legendre expansion for the metallic case.

Fig. 17
Fig. 17

Geometry of a checkerboard grating.

Fig. 18
Fig. 18

Reflection of order zero for the checkerboard grating shown in Fig. 17. Dielectric constants, 2 = 2.25, 1 = 5.76, out = 1. The heights are h1 = h2 = 10d. Note the high precision of the order-19 Legendre approximations for both polarization cases reflected by the exact coincidence of the resonances with that of the order-35 calculation.

Equations (78)

Equations on this page are rendered with MathJax. Learn more.

E ( x , y ) = E inc exp [ i ( α 0 x χ 0 y ) ] ,
α 0 = k 0 sin ( Θ ) , χ 0 = k 0 cos ( Θ ) .
k 0 = ( 2 π ) / λ .
[ x 2 + y 2 + k 0 2 ( x , y ) ] E ( x , y ) = 0 .
E p ± ( x , y ) = exp [ i ( α p x ± χ p y ) ] ,
α p = [ ( 2 π ) / d ] p + k 0 sin ( Θ ) , p = integer
χ p = χ p inc = inc k 0 2 α p 2 +
χ p = χ p out = out k 0 2 α p 2 +
α n = α p ( n ) , χ n = χ p ( n ) ,
E n ± ( x , y ) = ϕ n ( x ) exp ( ± i χ n y ) ,
ϕ n ( x ) = ( 1 / d ) exp ( i α n x )
E ( x , y ) = n = 1 [ exp ( i χ 1 inc y ) E inc δ n 1 + r n exp ( i χ n inc y ) ] ϕ n ( x )
E ( x , y ) = n = 1 t n exp ( i χ n out y ) ϕ n ( x )
R n = χ n inc χ 1 inc | r n E inc | 2
T n = χ n out χ 1 inc | t n E inc | 2
( x , y ) = ( j ) ( x ) .
E ( x , y ) = X ( x ) Y ( y )
X ( x ) Y ( y ) + X ( x ) Y ( y ) + k 0 2 ( j ) ( x ) X ( x ) Y ( y ) = 0 ,
X ( x ) X ( x ) + k 0 2 ( j ) ( x ) = Y ( y ) Y ( y ) .
X n ( x ) + k 0 2 ( j ) ( x ) X n ( x ) = χ n 2 X n ( x ) ,
X n ( x + d ) = exp ( i α 0 d ) X n ( x ) ,
E ( x , y ) = n = 1 N { a n + ( j ) exp [ + i χ n ( j ) y ] + a n ( j ) exp [ i χ n ( j ) y ] } X n ( j ) ( x ) ,
X n ( x ) = m = 0 M l p m ( l ) P m ( ξ ) , x l ( j ) x x l + 1 ( j ) ,
ξ = [ 2 x x l ( j ) x l + 1 ( j ) ] / [ x l + 1 ( j ) x l ( j ) ]
f ( x ) = l = 0 L a l P l ( x ) .
d 2 d x 2 f ( x ) = l = 0 L 2 b l P l ( x ) ,
b n = ( n + 1 2 ) p = n + 2 , n + 4 , . . . L ( p + n + 1 ) ( p n ) a p ,
P n ( ± 1 ) = ( ± 1 ) n , d d x P n ( x ) x = ± 1 = ( ± 1 ) n + 1 n ( n + 1 ) 2 ,
m = 0 M l p m ( l ) 1 = m = 0 m l + 1 p m ( l + 1 ) ( 1 ) m ,
α m = 0 M l p m ( l ) m ( m + 1 ) 2 = β m = 0 M l + 1 p m ( l + 1 ) ( 1 ) m + 1 m ( m + 1 ) 2 ,
α = 2 / { l ( j ) [ x l + 1 ( j ) x l ( j ) ] } ,
β = 2 / { l + 1 ( j ) [ x l + 2 ( j ) x l + 1 ( j ) ] } .
M = l = 1 L ( j ) ( M l 1 ) .
1 1 P n ( x ) exp ( i a x ) d x = i n 2 π a J n + ( 1 / 2 ) ( a ) .
M l = const . | l | [ x l + 1 ( j ) x l ( j ) ]
H ( x , y ) = n = 1 N { a n + ( j ) exp [ + i χ n ( j ) ( y y j ) ] + a n ( j ) exp [ i χ n ( j ) ( y y j ) ] } X n ( j ) ( x ) ,
1 ( j ) ( x ) y H ( x , y ) = 1 ( j ) ( x ) n = 1 N i χ n ( j ) × { a n + ( j ) exp [ + i χ n ( j ) ( y y j ) a n ( j ) exp [ i χ n ( j ) ( y y j ) ] } X n ( j ) ( x ) .
H ( x , y ) = n = 1 N { a n + ( k ) exp [ + i χ n ( k ) ( y y k ) ] + a n ( k ) exp [ i χ n ( k ) ( y y k ) ] } X n ( k ) ( x ) ,
1 ( k ) ( x ) y H ( x , y ) = 1 ( k ) ( x ) n = 1 N i χ n ( k ) × { a n + ( j ) exp [ + i χ n ( k ) ( y y k ) a n ( k ) exp [ i χ n ( k ) ( y y k ) ] } X n ( k ) ( x ) .
h j = y j + 1 y j
b n + = a n + ( j 1 ) exp [ + i χ n ( j 1 ) h j 1 ] ,
b n = a n ( j 1 ) exp [ i χ n ( j 1 ) h j 1 ] ,
X k ( j ) 1 ( j ) ( x ) X n ( j ) = δ k n .
a k + ( j ) + a k ( j ) = n = 1 N M k n ( j 1 ) ( b n + + b n ) ,
i χ k ( j ) [ a k + ( j ) a k ( j ) ] = n = 1 N i χ n ( j 1 ) N k n ( j 1 ) × ( b n + b n ) ,
M k n ( j 1 ) = X k ( j ) | 1 ( j ) ( x ) | X n ( j 1 ) ,
N k n ( j 1 ) = X k ( j ) | 1 ( j 1 ) ( x ) | X n ( j 1 ) .
a k + ( j ) = n = 1 N [ F k n ( j 1 ) b n + + G k n ( j 1 ) b n ] ,
a k ( j ) = n = 1 N [ G k n ( j 1 ) b n + + F k n ( j 1 ) b n ] .
F k n ( j 1 ) = 1 2 [ M k n ( j 1 ) + χ n ( j 1 ) χ k ( j ) N k n ( j 1 ) ] ,
G k n ( j 1 ) = 1 2 [ M k n ( j 1 ) χ n ( j 1 ) χ k ( j ) N k n ( j 1 ) ] .
a + ( j ) = F ( j 1 ) E ( j 1 ) a + ( j 1 ) + G ( j 1 ) 1 E ( j 1 ) a ( j 1 ) ,
a ( j ) = G ( j 1 ) E ( j 1 ) a + ( j 1 ) + F ( j 1 ) 1 E ( j 1 ) a ( j 1 ) ,
E k l ( j ) = exp [ i χ k ( j ) h j ] δ k l ,
d 2 d x 2 f ( x ) = 10 d d x f ( x ) + 11 f ( x )
f ( 1 ) = 1 , d d x f ( x ) x = 1 = 1 .
E 1 ( j ) a ( j ) = F 1 ( j ) a ( j + 1 ) F 1 ( j ) G ( j ) E ( j ) a + ( j ) .
a + ( j + 1 ) = [ F ( j ) G ( j ) F 1 ( j ) G ( j ) ] × E ( j ) a + ( j ) G ( j ) F 1 ( j ) a ( j + 1 ) ,
a ( j ) = E ( j ) F 1 ( j ) a ( j + 1 ) E ( j ) F 1 ( j ) G ( j ) E ( j ) a + ( j ) .
[ F ( 0 ) I 0 0 0 0 0 0 G ( 0 ) 0 I 0 0 0 0 0 0 I C 1 ( 1 ) C 2 ( 1 ) 0 0 0 0 0 0 C 3 ( 1 ) C 4 ( 1 ) I 0 0 0 0 0 0 I C 1 ( 2 ) C 2 ( 2 ) 0 0 0 0 0 0 C 3 ( 2 ) C 4 ( 2 ) I 0 0 0 0 0 0 I C 1 ( 3 ) 0 0 0 0 0 0 0 C 3 ( 3 ) I ] × [ a ( 0 ) a ( 1 ) a + ( 1 ) a ( 2 ) a + ( 2 ) a ( 3 ) a + ( 3 ) a + ( 4 ) ] = [ 0 0 0 0 0 0 C 2 ( 3 ) δ m , 1 C 4 ( 3 ) δ m , 1 ] .
C 1 = E F 1 G E ,
C 2 = E F 1 ,
C 3 = ( F G F 1 G ) E ,
C 4 = G F 1 ,
a ( L + 1 ) = E inc δ m 1 ,
f ( x ) = n = 0 a n T n ( x ) ;
d 2 d x 2 f ( x ) = n = 0 b n T n ( x ) ,
b n = 1 1 + δ n , 0 p = n + 2 , n + 4 , . . . p ( p 2 n 2 ) a p .
T n ( ± 1 ) = ( ± 1 ) n , d d x T n ( x ) x = ± 1 = ( ± 1 ) n + 1 n 2 .
1 1 T n ( x ) exp ( i a x ) d x 1 x 2 = i n π J n ( a ) .
X l | 1 ( x ) | X k = δ l k .
n X l | 1 ( x ) | ϕ n ϕ n X k = δ l k
n X l ϕ n ϕ n | 1 ( x ) | X k = δ l k .
n = 1 N X l adj | 1 ( x ) | ϕ n ϕ n X k = δ l k .
A ln = X l ( j ) adj | 1 ( j ) ( x ) | ϕ n ,
A n k ( j ) = ϕ n X k ( j ) .
M k l ( j 1 ) = X k ( j ) adj | 1 ( j ) ( x ) | X l ( j 1 ) = { [ A ( j ) ] 1 A ( j 1 ) } k l ,
N k l ( j 1 ) = X k ( j ) adj | 1 ( j 1 ) ( x ) | X l ( j 1 ) = { [ B ( j ) ] 1 B ( j 1 ) } k l .

Metrics