The corresponding Ewald summations are hybrid in the sense that the result is expressed in terms of two series:
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
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 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
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
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]]
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.
x | y | ReG | ImG | |
---|---|---|---|---|
D | 0.2 | 0.03 | 0.117120006144932 | -0.108131857633201 |
E | 0.117120006144932 | -0.108131857633206 | ||
YY | 0.117120006144941 | -0.108131857633206 | ||
NMcP | 0.117120006141860 | -0.108131857633197 | ||
D | 0.2 | 0.003 | 0.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 | ||
NMcP | 0.115881138135960 | -0.103450147416785 |
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].
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].
© Alexander Moroz, August 1, 2001 (last updated on September 25, 2010)