Exponentially convergent lattice sums of (quasi-) periodic Green's function of the Helmholtz and Laplace equations

Knowledge of the Green's function G of the scalar Helmholtz and Laplace equations in the so-called layer, or quasi-periodic, resp. quasiperiodic, case, i.e., when G is periodic in the number of dimensions that is smaller than the dimension of an embedding space, is a key to an efficient numerical solutions of boundary integral equations for potential flows in fluid mechanics, remote sensing of periodic surfaces, periodic gratings, and infinite arrays of resonators coupled to a waveguide, and, just to name a few, has applications in molecular dynamics and for the description of quasi-periodic arrays of point interactions in quantum mechanics.

An important class of problems which requires an efficient calculation of G in the Laplace case arises in many contexts of simulating systems of charged particles, such as crystal binding and lattice vibrations, Madelung constant, and in molecular dynamics in Monte Carlo simulations of particles interacting by long-range Coulomb forces. For instance, the periodic solution of the Laplace equation in two dimensions (2D) corresponds to the solution of particles interacting with the logarithmic interaction and has applications in simulations of 2D vortices in high-temperature superconductors.

For the analysis of diffraction of classical and quantum waves in various ab-initio first-principle multiple-scattering theories one then needs lattice sums, or the expansion coefficients D_L of G in the basis of regular [cylindrical in 2D and spherical in three dimension (3D)] waves. (See my recent review [24] for a list of further applications.)

However, even with the latest progress due to Yasumoto and Yoshitomi [1], it took 40 seconds to compute the Green's function G on SPARC workstation from lattice sums with 14 digits accuracy at a single point and frequency for the case of a one-dimensional (1D) periodicity in 2D. This was in striking contrast to the calculation of exponentially convergent lattice sums in the so-called bulk cases, i.e., when G is periodic in all space dimensions. The lattice sums for an infinite 2D lattice in 2D [2] and an infinite 3D lattice in 3D [3], which still seem to be largely unknown to some parts of optical and classical wave scattering communities, result from the so-called complete Ewald summation [4]. The corresponding Ewald summations are hybrid in the sense that the result is expressed in terms of two series, one involving summation over the spatial domain and the other over spectral domain [2,3]. The respective convergence times (on a PC with Pentium II processor) for a set of bulk 2D and 3D lattice sums with 6 digits accuracy are less than approx 0.03 second (for lmax=20) [5] and approx 0.8 second (for lmax=6) [6].

The complete Ewald summation is then a point of departure to derive the corresponding lattice sums, i.e., the expansion coefficients D_L of G in the basis of regular [cylindrical in 2D and spherical in three dimension (3D)] waves. Using an extension of the Ewald summation procedure, exponentially convergent lattice sums were also derived in the quasiperiodic, case. For 2D periodicity in 3D, exponentially convergent expressions for G and its lattice sums were derived in the late sixties of the last century by Kambe [5-6]. Since then they became an integral and indispensable part of numerical programs, based on the layer Korringa-Kohn-Rostoker (LKKR) method [7-15] within low-energy electron-diffraction (LEED) theory [11,16], which describe scattering and diffraction of quantum and classical waves off (a finite stack of) 2D lattice(s) in 3D. The Kambe results are so inherently built in into the LEED/LKKR theory that nowday they are hardly ever cited when some results are obtained by means of the LEED/LKKR theory. The LKKR and LEED based numerical programs have been successfully tested, time and again and for various physical problems, since the early seventies of the last century. The last decade they have also been incorporated into acoustic, elastic, and electromagnetic variants of the LKKR theories [12-15]. The convergence of Kambe's expressions has been recently studied in Ref. [21] and it was shown to be superior to some recently derived alternative expressions of the lattice sums.

Yet, an extension (or reduction if you wish) of the Kambe analysis to one dimension lower, i.e., to the layer case of 1D periodicity in 2D has not been performed till very recently [17]. An extension of the Kambe analysis to the layer case of 1D periodicity in 3D has been only recently presented in my article on Quasi-periodic Green's functions of the Helmholtz and Laplace equations [math-ph/0602021]. Yes, there have been several attempts to calculate the lattice sums [1,18]. However, none of them succeeded in providing an exponentially convergent analytic representation of lattice sums either for 1D periodicity in 2D or for 1D periodicity in 3D. Moreover, the new representations [see e.g. Eqs. (7)-(9) of my article Exponentially convergent lattice sums, Opt. Lett. 26, 1119 (2001) [pdf]]

• (i) can be implemented numerically more simply;
• (ii) converge roughly 200 times faster than the previous best representation [1];
• (iii) can be directly incorporated (see Ref. [19]) into powerful LKKR [7-15] and LEED [11,16] techniques developed within the context of diffraction and scattering of electron waves, which, such as illustrated in Refs. [5,6], can also be applied to a variety of classical, including electromagnetic, waves.

Numerical evaluation: As a test of numerical efficiency, the case of a plane wave incident at an angle pi/8 upon 1D grating oriented along the x-axis was investigated. In this case, for wavelength lambda/v_0=0.23 and fixed x/v_0=0.2 (v_0 is the length of a period (primitive lattice cell) along the x-axis), convergence results of several other methods for G have been presented in the literature for different values of y/v_0 (see Table IV of Ref. [1] and Table III of Ref. [18]). Results obtained by my method and their comparison with those obtained by previous methods are summarized in Table I. The cutoff value lmax on the summation over l in Eq. (4) of Opt. Lett. 26, 1119 (2001) depends on the product s=sigma*Rmax, where sigma=2*pi*v_0/lambda and Rmax is the largest of |R| for which G is going to be calculated. Since Bessel functions in Eq. (4) of Opt. Lett. 26, 1119 (2001) become rapidly decaying for l>s, a value of lmax slightly larger than s is usually enough. In the current case, to achieve a comparable accuracy to that reported in Table IV of Ref. [1] and Table III of Ref. [18], the cutoff lmax was taken to be lmax=50, i.e., significantly larger than s. The cutoff lmax=50 was also used in Refs. [1,18]. Bessel functions were generated with the help of double-precision routines bessjy, beschb, and chebev of Numerical Recipes in Fortran 77 (see Sec. 6.7. there). Initial recurrence values of definite integrals in Eqs. (8) and (10b) of Opt. Lett. 26, 1119 (2001) were performed using routines qrombt, polint, and trapzd of Numerical Recipes in Fortran 77. The routine qromb was used with the parameters:

An initial recurrence value of the incomplete gamma function in Eq. (8) of Opt. Lett. 26, 1119 (2001) for n=0 and for a real second argument, which is equivalent to the error function erfc [see Eq. (10a) of Opt. Lett. 26, 1119 (2001)], was calculated using routines erfc, gammp, gcf, gser, gammln, and gammq of the same book. Convergence parameter EPS in the routines gcf and gser was set to EPS=1.E-15.

The Ewald parameter eta in Eqs. (7)-(9) of Opt. Lett. 26, 1119 (2001) can often be varied several order of magnitude without affecting the result in a wide frequency window. However, for some singular values of eta, one can enter a numerically instable region: the D_l^{(1)} and D_l^{(2)} contributions have opposite sign and similar magnitude, which is several orders larger than the resultant D_l [Eqs. (6) of Opt. Lett. 26, 1119 (2001)]. This instability can be easily remedied by choice of some other value of eta, or, one can make eta depend on sigma and l and prevent numerical instability completely [20]. Of the cases tested, the simplest case of a constant eta was chosen, as was the case in Refs. [5,6]; here the numerical instability limited limited eta to the interval (0, 0.2). The numerical code OLA for a 1D periodicity in 2D is available here.

Convergence time: Having achieved comparable precision to that of Yasumoto and Yoshitomi [1] and Nicorovici and McPhedran [18] (see Table I), respective computational times could have been compared for lmax=50. Only the convergence time T1 (see Tab. III of Ref. [18]), i.e., the time needed to calculate G at a single point, was investigated. A reliable comparison of convergence times to calculate G at several points with those reported in Table IV of Ref. [1] and Table III of Ref. [18] would require the use of identical routines to calculate the Bessel functions. In our case, the computational time to reproduce a value of G within 8.E-15 of that obtained by a direct summation turns out to be approx 0.2 second (F77 program was compiled without any optimalization), in line with the respective approx 0.03 and 0.8 second for the convergence time of a set of bulk 2D [5] and 3D lattice sums [6] with 6 digits accuracy. This should be compared to 40 seconds of Yasumoto and Yoshitomi [1] and to 1232 seconds of Nicorovici and McPhedran [18]. The direct summation of the dual representation (spectral domain form) of G (the entries in rows labelled by D in Table I) took 4 seconds for y=0.03, 36 seconds for y=0.003, and 366 seconds for y=0.0003 (see Tab. III of Ref. [18]). Therefore, even for a single point, the expressions (7)-(9) of Opt. Lett. 26, 1119 (2001) are also superior to variants of the spectral domain representation of G. (Only for a very large y, (a variant of) the spectral domain representation will have the same speed of convergence as the method based on the Eqs. (7)-(9) of Opt. Lett. 26, 1119 (2001), since its convergence becomes also exponential.) If, given sigma and k_parallel, G is needed at several points, this advantage is enhanced further. Indeed, let as above be T1 the time needed to calculate G at a single point. Then, after initial calculation of the lattice sums D_l up to a cutoff value lmax, any further evaluation of G for a given sigma and k_parallel only requires a straightforward evaluation of regular Bessel functions and performing the sum in Eq. (4) of Opt. Lett. 26, 1119 (2001) with the same set of the D_l's. The latter cost only a tiny fraction (cca 1/100) of T1. On the other hand, using variants of the spectral domain representation of G, the time needed to calculate G at N points would be N*T1. Therefore, compared to variants of the spectral domain representation of G, expressions (7)-(9) of Opt. Lett. 26, 1119 (2001) provide an unparallel speed of convergence for G.

In analogy to the bulk case [22], the results for lattice sums in Opt. Lett. 26, 1119 (2001) can easily be generalized to any complex (non-Bravais) lattice, as it was accomplished by Kambe [23] for the case of a 2D periodicity in 3D. Numerical F77 codes implementing Kambe's formulas for the case of a 2D periodicity in 3D can be found in a book by Pendry [11], or, can be downloaded from Comp. Phys. Commun. [25,26]. Together with those for a 1D periodicity in 2D they are also available upon request from the author.

Following the link rta1in2k.exe you can download a Windows executable, which employes the lattice sums up to the angular-momentum cutoff lmax=15, and which calculates reflection, transmission, and absorption off a photonic crystal of a finite width formed by stacking of 25 layers along the X-direction of a triagonal lattice of coated cylindrical rods. For more results see Ref. [19].

References

1. K. Yasumoto and K. Yoshitomi, Efficient calculation of lattices sums for free-space periodic Green's function, IEEE Trans. Antennas Propag. 47(6), 1050-1055 (1999).
2. A. M. Ozorio de Almeida, Real-space method for interpreting micrographs in cross-grating orientations. I. Exact wave-mechanical formulation, Acta Cryst. A 31, 435-442 (1975) (the resulting formulae for k=0 contain a minor misprint).
   [For a generic for k see J. S. Faulkner, Multiple scattering calculation in two dimensions, Phys. Rev. B 38, 1686-1694 (1988) and Eqs. (27) therein,
   and Appendix of the article by K. M. Leung and Y. Qiu, Multiple-scattering calculation of the two-dimensional photonic band structure, Phys. Rev. B 48, 7767–7771 (1993).]
3. F. S. Ham and B. Segall, Energy bands in periodic lattices - Green's function method, Phys. Rev. 124, 1786-1796 (1961).
4. P. P. Ewald, Die Berechnung optischer und elektrostatistischer Gitterpotentiale, Ann. Phys. (Leipzig) 64, 253-287 (1921).
5. H. van der Lem and A. Moroz, Towards two-dimensional complete photonic-bandgap structures below infrared wavelengths,
   J. Opt. A: Pure Appl. Opt. 2, 395-399 (2000);
   

   H. van der Lem, A. Moroz, and A. Tip, Band structure of absorptive two-dimensional photonic crystals,
   J. Opt. Soc. Am. B 20, 1334-1341 (2003) [pdf].

6. A. Moroz and C. Sommers, Photonic band gaps of three-dimensional face-centered cubic lattices,
   J. Phys.: Condens. Matter 11, 997-1008 (1999) [physics/9807057] [story behind this article].
7. K. Kambe, Theory of Electron Diffraction by Crystals. Green's function and integral equation, Z. Naturforschg. 22a, 422-431 (1967). See also E. G. McRae, Multiple-scattering treatment of LEED intensities, J. Chem. Phys. 45, 3258-3276 (1966). .
8. K. Kambe, Theory of Low-Energy Electron Diffraction. I. Application of the Cellular method to monoatomic layers, Z. Naturforschg. 22a, 322-330 (1967).
9. P. J. Jennings and E. G. McRae, Electron diffraction at crystal surfaces IV. Computation of LEED intensities for ``muffin-tin" models with application to tungsteen (001), Surf. Sci. 23, 363-388 (1970).
10. D. W. Jepsen, P. M. Marcus, and F. Jona, Accurate calculation of the low-energy electron-diffraction spectra of Al by the layer-Korringa-Kohn-Rostoker method, Phys. Rev. Lett. 26, 1365-1368 (1971).
11. J. B. Pendry, Low Energy Electron Diffraction (Academic Press, London, 1974).
12. K. Ohtaka, Scattering theory of low-energy photon diffraction, J. Phys.: Solid State, Vol. 13, 667-680 (1980).
13. A. Modinos, Scattering of electromagnetic waves by a plane of spheres-formalism, Physica 141, 575-588 (1987).
14. N. Stefanou, V. Karathanos, and A. Modinos, Scattering of electromagnetic waves by periodic structures, J. Phys.: Condens. Matter 4, 7389-7400 (1992).
15. I. E. Psarobas, N. Stefanou, and A. Modinos, Scattering of elastic waves by periodic arrays of spherical bodies, Phys. Rev. B 62, 278 (2000).
16. E. G. McRae, Electron diffraction at crystal surfaces I. Generalization of Darwin's dynamical theory, Surf. Sci. 11, 479-491 (1968).
17. A. Moroz, Exponentially convergent lattice sums, Opt. Lett. 26, 1119-1121 (2001) [pdf] (accompanying code OLA).
    [See also Appendices A and B of the article by K. Ohtaka, T. Ueta, and K. Amemiya, Calculation of photonic bands using vector cylindrical waves and reflectivity of light for an array of dielectric rods, Phys. Rev. B 57, 2550–2568 (1998).]
18. N. A. Nicorovici and R. C. McPhedran, Lattice sums for off-axis electromagnetic scattering by gratings, Phys. Rev. E 50, 3143-3160 (1994).
19. V. Poborchii, T. Tada, T. Kanayama, and A. Moroz, Silver-coated silicon-pillar photonic crystals: enhancement of a photonic band gap,
    Appl. Phys. Lett. 82, 508-510 (2003) [pdf];
    A. Moroz, On layer Korringa-Kohn-Rostoker method in two dimensions, in preparation.
20. M. V. Berry, Quantizing a classically ergodic system: Sinai's billiard and the KKR method, Ann. Phys. (NY) 131, 163-216 (1981).
21. A. Moroz, On the computation of the free-space doubly-periodic Green's function of the three-dimensional Helmholtz equation,
    J. Electromagn. Waves Appl. 16, 457-465 (2002).
22. B. Segall, Calculation of the band structure of ``complex" crystals, Phys. Rev. 105, 108-115 (1957).
23. K. Kambe, Theory of Low-Energy Electron Diffraction. II. Cellular method for complex monolayers and multilayers, Z. Naturforschg. 23a, 1280-1294 (1968).
24. A. Moroz, Quasi-periodic Green's functions of the Helmholtz and Laplace equations,
    J. Phys. A: Math. Gen. 39, 11247-11282 (2006). [erratum] [math-ph/0602021]
25. J. M. McLaren, S. Crampin, D. D. Vvedensky, R. C. Albers, and J. B. Pendry, Layer KKR electronic structure code for bulk and interface geometries, Comp. Phys. Commun. 60, 365-389 (1990).
26. V. Yannopapas, N. Stefanou, and A. Modinos, Heterostructures of photonic crystals: frequency bands and transmission coefficients, Comp. Phys. Commun. 113, 49-77 (1998).