Efficient calculation of lattice sums and Green's functions * Ewald summation * Exponentially convergent lattice sums * Free-space periodic and subperiodic, quasi-periodic, or quasiperiodic, Green's function of the Helmholtz and Laplace equations * Schlömilch-type series * Schloemilch series

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

Ewald summation

The knowledge of a periodic Green's function G of scalar Helmholtz and Laplace equations is of key importance to an efficient numerical solution of numerous important problems in classical physics and quantum mechanics. An efficient calculation of the periodic Green's function G of the scalar Helmholtz equation can be accomplished by means of the complete Ewald summation as shown

The corresponding Ewald summations are hybrid in the sense that the result is expressed in terms of two series:

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].

A particular class of problems requires the knowledge of Green's function G in the so-called layer, subperiodic, 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. Such a Green's function G is 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. (See my recent review [24] for an overview of applications and list of references.)

However a subperiodic Green's function G is much more difficult to obtain than the periodic Green's function. Whereas, in the latter case, the Ewald summation requires nothing than an eigenfunction expansion of the Green's function with a subsequent application of a generalized Jacobi identity [2-4,20,24], the Ewald summation is much more involved in the former case [7,8,17,24]. In the subperiodic case of 2D periodicity in 3D the Ewald summation has been performed for the very first time by

both in connection with the applications involving the so called low-energy electron-diffraction (LEED) theory [11,16].

Unfortunately, these developments within the LEED theory went for a very long time unnoticed by the various classical wave communities, such as electromagnetic community etc. By that I mean the time that defies any imagination, i.e. sometime longer than 40 years!!! For instance, the electromagnetic community still continues to celebrate the article by

as the article that for the first time provided the Ewald summation for the subperiodic case of 2D periodicity in 3D. Yet, at that time, already the more advanced subject of exponentially convergent lattice sums had long time been described in monographs such as Low Energy Electron Diffraction: The Theory and Its Application to Determination of Surface Structure by J. B. Pendry that was published as early as in 1974! Nevertheless, ``searching" for the ``solution" of this already solved problem would continue to occupy many researches from various classical wave communities well into the 21th century. I wonder how much money of various grant agencies has been burn out and how many people have got their tenures in the process of rediscovering the known solutions? As a rule of thumb, if you find a paper that deals with either Ewald summation or lattice sums in a subperiodic case but does not quote the work of Kambe, you can be almost certain that the author(s) did not make his(their) homework.

It appears that two extraterrestrial civilizations would find easier and earlier a common language to communicate than the classical wave community will attempt to understand the results of quantum wave community that uses the very same Helmholtz equation to describe wave propagation!!!

Lattice sums

For an efficient analysis of diffraction of classical and quantum waves in various ab-initio first-principle multiple-scattering theories the so-called lattice sums D_L are required. The latter are defined as the expansion coefficients of G in the basis of regular [cylindrical in 2D and spherical in 3D] waves. With the help of the expansion coefficients D_L it is required to perform only a single summation over the lattice in order to get the values of G in all spatial points. This distingushing feature makes the use of D_L superior to other techniques which, as a rule, require repeated lattice summation for each spatial point in order to determine G in that point. Once the scalar lattice sums are determined, the lattice sums for the remaining vectorial and tensorial problems can be obtained by multiplying the scalar lattice sums with suitable geometrical coefficients. (See my recent review [24] for further details and references.)

Lattice sums for a periodic Green's function G, which still seem to be largely unknown to some parts of optical and classical wave scattering communities, can be rather easily obtained:

The point of departure and a prerequisite step in order 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 3D) waves, is the complete Ewald summation. However, mirroring the Ewald summation, to obtain lattice sums for a subperiodic Green's function G is much more involved. Using an extension of the Ewald summation procedure, exponentially convergent expressions for the lattice sums D_L of G in the subperiodic case of 2D periodicity in 3D were for the very first time obtained after elaborate and tedious calculations of

in the late sixties of the last century. Like in the case of the lattice sums for bulk Green's functions, the resulting Kambe's expressions do not contain any Bessel function: they all disappear in a limiting step.

Since then, Kambe's result had became an integral and indispensable part of numerical programs based on the layer Korringa-Kohn-Rostoker (LKKR) method [7-15] within the LEED theory [11,16], which describe scattering and diffraction of quantum and classical waves off (a finite stack of) 2D lattice(s) in 3D. See in this regard the monograph by

that in an Appendix provides a full listing of Fortrane (F77) LEED code. The routines of Pendry's LEED code can be used even by a pragmatic engineer as black box routines. The Pendry's LEED code comprises among others also, up to minor amendmends, the program lines that I have arranged in a routine dlsumf2in3.f for the calculation of Kambe's lattice sums. My code gff2in3.f can be then used as a main code. It generates the lattice points, calls the routines that calculate spherical harmonics and spherical Bessel functions, and helps you to obtain the free-space Green's function given the total lattice sum D_L generated by dlsumf2in3.f.

The total lattice sum D_L is conventionally written as a sum of three partial contributions

D_L = D_L^1 + D_L^2 + D_L^3

If you look at my J. Phys. A: Math. Gen. 39, 11247-11282 (2006), the D_L^1 term for different subperiodic cases is given by Eqs. 83,85,86. The only special function is the incomplete gamma function which is easily calculated numerically on using recurrences given by Eqs. 88-89. There is nothing mysterious or difficult in calculating incomplete gamma function compared to e.g. calculation of the Bessel functions.

The D_L^2 term for different subperiodic cases is given by Eqs. 102-103. There is no special function involved. An integral therein is calculated easily numerically again on using recurrence given by Eq. 106.

The D_L^3 term for different subperiodic cases is given by Eqs. 118-119. There is no special function involved and the sum is calculated numerically in a straightforward way.

The Kambe results have been so inherently built in into the LEED/LKKR theory that nowday they are hardly ever cited when some new results obtained by means of the LEED/LKKR theory are reported. 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.

Kambe's effort had really payed off. Without the use of Kambe's techniques, 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.

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]. The latter article provides a unified treatment of all the subperiodic cases (i.e. 1D and 2D periodicity in 3D, together with 1D periodicity in 2D) within the spirit of Kambe's approach and shows how to derive all the results for the separate lattice sums in a single go.

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]]

In addition to the applications listed above (see also my recent review [24]), the results presented therein make it in principle possible efficient investigations (a) of the diffraction of light by (photonic crystal) gratings and (b) of subtle phenomena involving light-matter interaction of light in the presence of a single or a finite stack of 1D gratings [19].

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:

PARAMETER (EPS=1.d-15, JMAX=30, JMAXP=JMAX+1, K=10, KM=K-1)

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.

Table I. Free-space periodic Green's function G for off-axis incidence at an angle theta=pi/8 upon a one-dimensional lattice oriented along the x-axis in two-dimensions with lambda/v_0=0.23. Here v_0 is the length of a period (primitive lattice cell) along the x-axis. In rows labeled by D, values of G obtained by a direct summation of its dual representation (spectral domain form) (taken from Tab. III of Ref. [18]). These data are compared against those in rows labeled by E, obtained by the (complete) Ewald summation [Eqs. (7)-(9) of Opt. Lett. 26, 1119 (2001)] with the Ewald parameter eta=0.011. Data in the respective rows labelled by YY and NMcP are those obtained by Yasumoto and Yoshitomi [1] and Nicorovici and McPhedran [18] methods.

x y ReG ImG
E 0.117120006144932-0.108131857633206
NMcP 0.117120006141860-0.108131857633197
D 0.20.0030.115891895634567-0.103497063599642
E 0.115891895634565-0.103497063599651
YY 0.115891895634577-0.103497063599646
NMcP 0.115891895630095-0.103497063599643
D 0.2 0.0003 0.115881138140449-0.103450147416784
E 0.115881138140448-0.103450147416794
YY 0.115881138140457-0.103450147416788

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 that 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].


  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 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: The Theory and Its Application to Determination of Surface Structure (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).

[ Main page | My personal webpage | Full list of my scientific publications | List of my publications in Scitation | My most frequently cited articles]
[ Photonic Crystals Headlines 25 Years Old, at Least ... | Selected Links on Photonics, Photonic Crystals, Numerical Codes, Free Software
[ Negative refractive index material headlines long before Veselago work ]

© Alexander Moroz, August 1, 2001 (last updated on September 25, 2010)

eXTReMe Tracker