Choosing Splitting Parameters and Summation Limits in the Numerical Evaluation of 1-D and 2-D Periodic Green's functions Using the Ewald Method

[ 1 ] Accurate and efficient computation of periodic free-space Green’s functions using the Ewald method is considered for three cases: a 1-D array of line sources, a 1-D array of point sources, and a 2-D array of point sources. A limitation on the numerical accuracy when using the ‘‘optimum’’ E parameter (which gives optimum asymptotic convergence) at high frequency is discussed. A ‘‘best’’ E parameter is then derived to overcome these limitations. This choice allows for the fastest convergence while maintaining a specific level of accuracy (loss of significant figures) in the final result. Formulas for the number of terms needed for convergence are also derived for both the spectral and the spatial series that appear in the Ewald method, and these are found to be accurate in almost all cases.


Introduction
[2] In applying numerical full wave methods like the Method of Moments (MoM) or Boundary Integral Equations (BIE) to periodic structures involving conducting or dielectric electromagnetic scatterers, fast and accurate means for evaluating the free-space periodic Green's function (FSPGF) are often needed. This type of Green's function arises in a wide variety of applications, ranging from microwaves to optics, to the study of metamaterials and nanostructures. The Ewald method is one of the fastest methods for calculating the FSPGF. In the Ewald method, the FSPGF is expressed as the sum of a ''modified spectral'' and a ''modified spatial'' series. The terms of both series possess Gaussian decay, leading to an overall series representation that exhibits a very rapid convergence rate. The convergence rate is optimum when the ''optimum'' value of the Ewald splitting parameter is used [Jordan et al., 1986], denoted here as E opt . (For some applications, such as when using a periodic Green's function in a MoM solution with full-domain basis functions, one may not wish to have a balanced convergence between the two series, as explained in Mathis and Peterson [1996] and Mathis and Peterson [1998]. However, when using subdomain basis functions and performing the integrations over the basis and testing functions in the spatial domain, the objective is to minimize the computation time of the periodic Green's function, and this is accomplished by using E opt . In this case the singular integrals involved in evaluating the matrix elements can be handled by specially designed numerical quadrature rules [Khayat and Wilton, 2005].) [3] However, the numerical accuracy of the Ewald method degrades very quickly [Kustepeli and Martin, 2000] with increasing frequency (i.e., when the periodicity becomes large relative to a wavelength). This is due to a catastrophic loss of significant figures in combining the contributions of the two series, wherein the leading terms (and to a lesser extent, other nearby terms) in each series become very large but nearly equal and opposite in sign.
[4] The method proposed and studied here limits the size of the largest terms in the series relative to that of the total Green's function by modifying the value of the splitting parameter E to avoid undue loss of accuracy. Increasing the E parameter limits the size of the largest terms in both series at the expense of decreasing the convergence rate. Hence, there is a tradeoff between the size of the largest term allowed, which determines the number of significant figures lost, and the series convergence rate. A value E L of the Ewald splitting parameter is then obtained based on the number of significant figures L that may be lost. This ''best'' value, E L , then yields the fastest convergence of the Ewald series while limiting the loss of significant figures to the user-defined value L.
[5] A preliminary and intuitive analysis for the ''best'' choice of the Ewald splitting parameter was performed in Capolino et al. [2005] for the case of 1-D array of line sources, and in Capolino et al. [2007] for the case of a 1-D array of point sources. The case of a 2D-array of point sources has been treated in detail in Oroskar et al. [2006]. Here we extend the analysis of the last paper to the two other cases, and provide a unified formalism for the choice of the Ewald splitting parameter and the summation limits that is valid for all three cases.

Spectral, Spatial, and Ewald Green's Function Representations
[6] Three different FSPGF cases are considered: a 1-D array of line sources with interelement period d, a 2-D array of point sources on a general skewed lattice, and a 1-D array of point sources with interelement spacing d. The geometries are depicted along with relevant coordinate systems and geometrical definitions in Figures 1, 2, and 3, respectively. An interelement phase shift along the array direction(s) is assumed. Both spatial and spectral representations of each Green's function exist in the general form where the terms of the spatial representations are whereas those of the spectral representations arẽ 2jk zp e Àjk zp Dz ; 1 À D array of line sources 1 À D array of point sources where Dz |z À z 0 | and an e jwt time dependence is assumed and suppressed. In the above, (x, y, z) and (x 0 , y 0 , z 0 ) denote the observation and reference-element source points, respectively. For the 1-D cases, the periodicity is along x with a period of d (see Figures 1 and 2), and For the 2-D array with lattice vectors s 1 , s 2 (see Figure 3), the parameters are defined as The transverse phasing wave vector k t00 =xk sin q 0 cos f 0 + yk sin q 0 sin f 0 defines the interelement phasing for the 2-D array in terms of the propagation angles q 0 , f 0 of the (0, 0) Floquet mode. For 1-D arrays, this quantity becomes the scalar phasing k x0 = k cos q 0 , where q 0 is measured with respect to the x-axis. Physically the FSPGF is the timeharmonic scalar potential produced by the associated array of phased sources.
[7] For a lossless medium the square root for the wave numbers k rp and k zpq is that which gives a positive real number or a negative imaginary number, corresponding to waves that propagate outward or decay away from the sources, respectively. For a lossy medium the wave numbers are complex and it suffices to require the imaginary part of the wave numbers to be negative. Henceforth, it will be assumed that the medium is lossless, but the formulas can be extended to the lossy case.
[8] When employing the Ewald method for the evaluation of the FSPGF, the Green's function is expressed as a sum of two series [Ewald, 1921] of the form   q. The terms of the modified spectral representations arẽ For the 1-D array of point sources, the exponential integral function that appears may have a negative argument, depending on the frequency. In this case the argument is interpreted as being infinitesimally above the branch cut of the exponential integral function on the negative real axis, corresponding to an infinitesimal amount of loss.

Splitting Parameter (E)
[9] In equations (5) and (6) the spatial and spectral series both involve a ''splitting'' parameter E. The ''optimum value'' E opt for the splitting parameter [Jordan et al., 1986] balances the asymptotic rate of convergence of the spatial and spectral series, and consequently minimizes the overall number of terms needed to calculate the total FSPGF. The optimum value is found to be However, numerical difficulties are encountered when the lattice separations (periods) become large relative to a wavelength. This was first discovered in Kustepeli and Martin [2000], and subsequently also discussed in Capolino et al. [2005], Oroskar et al. [2006], and Capolino et al. [2007]. This happens because, for large arguments, both the complementary error function and the exponential integrals contribute terms of the form exp[(k/(2E)) 2 ] that produce extremely large initial terms (and to a lesser extent, large nearby terms) in both the spatial and the spectral series. The resulting series then each converge to very large, but nearly equal and oppositely signed values, resulting in a total sum of moderate value but with a catastrophic loss of significant figures when the two series are combined. Exponential overflow is another potential concern as well, due to the large initial values in each series.
[10] To circumvent the problem, it is desirable to limit the size of the largest terms of each series by choosing an E value larger than the ''optimum'' value, which reduces the maximum values of both the complementary error function and the exponential integrals. As a result, one avoids loss of accuracy in adding the two series, and a more accurate result for the total Green's function is obtained [Kustepeli, and Martin, 2000] at the expense of slower convergence.
CELEPCIKAY ET AL.: NUMERICAL EVALUATION OF EWALD METHOD 4 of 11

RS6S01
[11] In the following sections, a recipe for finding the best choice for E, called E L , that achieves the fastest convergence under the constraint of limiting the loss of significant figures to L digits, is obtained for a general source and observation point, for all three cases.

Choice of the Splitting Parameter
[12] The strategy is to limit the size of the largest terms relative to the value of the total Green's function, with the largest terms in each series being the initial (0) or (0, 0) terms in 1-D or 2-D, respectively. The value of E = E L is obtained by enforcing the following conditions: where G est is a closed-form estimate of the FSPGF, and the integer parameter L indicates (roughly) the number of significant figures one is willing to sacrifice in the calculation. For a strict bound, the factor a in (8) should be chosen as 1/2 to ensure that each of the two initial terms in (8) (one from the spatial series and one from the spectral series) contributes no more than half the total error limit. However, it is found in each of the three cases that one of the initial terms (either the spatial one or the spectral one, depending on the case) is significantly larger than the other, and a factor of a = 1 therefore becomes more appropriate, and this is adopted here. It suffices to use a rough estimate of the Green's function on the right-hand side of (8), which may be obtained by examining the most nearly singular terms from both the spatial and the spectral series. (This is an improvement over the previous derivations, as in Oroskar et al. [2006], where only the spatial term was used.) This yields the following magnitude estimate for the overall Green's function: In (9), the indices p, q, m, n that produce the smallest values for R m , R mn , k rp , k zp , and k zpq must be determined; although this may be done analytically, it is also very easy to find these values from a simple numerical search.
[13] To proceed with the derivation of E L , we replace the terms on the left hand side of (8) by their highfrequency asymptotic estimates. The complementary error function terms are estimated by using the asymptotic relation valid for large arguments. In addition, in (5) and (6) series of exponential integral functions appear. The simplest way to asymptotically evaluate these series is to re-cast them back into their integral forms. The integral forms are equation (15) in Capolino et al. [2007] and equation (11) in Capolino et al. [2005], which are reproduced below: The integrals that appear in the above identities may be asymptotically evaluated for k!1 via integration by parts [Felsen and Marcuvitz, 1994;Bleistein and Handelsman, 1986] (see case 1a and case 2a of Appendix A for further details), yielding Applying (9) - (14) in (8), as detailed in Oroskar et al.
[2006] yields in each case (spectral and spatial) a transcendental equation of the form where the constant K and the function F in (15) are defined in Table 1 (in each case the function F contains a factor c 1 that is also shown). Solving the quadratic form for z on the left-hand side of (15) puts the equation into the following form that may be efficiently solved iteratively, due to the slow variation of the ln function: where F i = F(z i ).
[14] The solution of (15) yields the parameter z, which (see Table 1) is inversely proportional to the desired Ewald parameter E. Two E values, E spectral and E spatial , result from this procedure. To properly bound both series, E should then be chosen as The resulting value is the smallest value of E that ensures that the largest term in both the spectral and the spatial series (see equation (8)) is limited in magnitude to avoid losing more than L significant figures when the two series are combined. At low or moderate frequency, where the initial terms of the spectral and spatial series are not large, E L = E opt since E opt will be the largest of the three terms in (17).

Number of Terms Needed for Convergence
[15] Having determined the ''best'' value of the E parameter E L as a function of frequency, the next goal is to determine how many terms should be summed in each series (modified spatial and modified spectral) in (4) to achieve convergence. (One could, of course, check convergence as the series are summed, but we eventually want to select between the Ewald and alternative methods for computing the FSPGF that may be more efficient in some situations; for that an a priori estimate of the number of series terms and their relative computational cost is required.) Recall that a given value of L, the number of significant figures sacrificed in the calculation, has been assumed. For the resulting value of E = E L we must then determine how many terms in each series are needed to guarantee convergence of the Green's function to S significant figures. A method is developed here to calculate the series index limits ±P, ±Q, ±M, and ±N for the series indices p, q, m, and n, respectively (Q and N only apply for the 2-D geometry). If the accuracy of the arithmetic is limited to T significant figures due to the machine precision (or, more likely in practice, limited by the accuracy of the complementary error function and exponential integrals), the value of S specified should be limited to S < T À L in order to avoid unnecessary computation.
[16] Owing to the Gaussian convergence, a rough estimate of the truncation error for both the spectral and spatial series is obtained by using the sum of the largest of the first neglected terms along each principal sum index in the series. Limiting the convergence error in each series to one half of the total, we thus require that The error in stopping the summations is approximated in the above equations as the sum of the two (four) values that give the largest contributions to the summed values for 1-D (2-D) arrays. The factor b is introduced to allow for an empirical adjustment of the summation limits. Using b = 1 corresponds to a strict error bound based on the assumed asymptotic approximations, and therefore usually represents a worst-case error bound. However, factors of b = 2 and b = 4 have been found to work well for the 1-D and 2-D cases, respectively. If we assume that the contributions of the positively indexed terms on the LHS of inequalities (18) and (19) are dominant, the choice of the limits amounts to choosing M, N, P, Q such that (The dominance of the positively indexed terms may be assumed without loss of generality, by placing absolute values on the spatial displacements and phasing wave numbers, as seen in (22) -(25).) Using the asymptotic estimates of Table A1, case 1b and case 2b, which assume large values of M, N and P, Q for the terms in (20) above, yields a transcendental equation of the form where D = x 2 + c 2 2 , and the constants c 2 and W in (21) are defined in Table 2.
[17] Equation (21) must be solved numerically for x, which, as seen in Table 2, is the only term involving the indices. Equation (21) may be efficiently solved iteratively for the parameter D as outlined in Oroskar et al. [2006].
[19] The results for the case of a 1-D line source array are For the 1-D point source array the results are For the 2-D point source array the result is (restricting the result to the case of the rectangular lattice s 1 = ax, s 2 = bŷ for simplicity): M ¼ Int In these final results, noninteger solutions of (21) for the summation limits should be rounded up to the next larger integer to obtain M + 1, etc., or equivalently, rounded down to obtain M, etc., as assumed in (22) -(25), in which the Int function truncates to the next lower integer.
[20] One note regarding (22) -(25) should be made in connection with the square roots. Depending on the geometry of the problem and the specified convergence accuracy, it may occur that the argument of one of the square roots is negative, yielding a complex value for the summation limit. This occurs because the asymptotic approximation used to estimate term magnitudes becomes invalid when the series actually needs only a few terms to converge, corresponding to a summation limit of zero or one. The problem is circumvented by always using a summation limit that is at least unity.

Results
[21] In this section results are presented for the three cases: 1-D line sources, 1-D point sources, and 2-D point sources. For all results, free-space conditions are assumed (k = k 0 and l = l 0 ). For the 1-D cases, d = 0.5 m. The 2-D results are shown for a square lattice (a = b = 0.5 m). As a result, the optimum splitting parameter is E opt = 3.5449 for all three cases. A zero progressive phase shift is assumed (all source elements are in phase). The reference source element is located at the origin and the observation point is located at (0, 0, Dz) with the vertical distance from the source plane set at Dz = 0.05 m. For brevity's sake, only the magnitude of the Green's function terms are shown in the results.
[22] For each case, three tables are shown. For the 1-D array of line sources, Tables 3a -3c are shown; for the 1-D array of point sources, Tables 4a -4c are shown; and for the 2-D array of point sources, Tables 5a -5c are shown. The first tables in each case (part (a)) illustrate that the values ofG 0,0 E and G 0,0 E become enormous at high frequency. For these tables the number of lost significant digits is set to L = 3. (The calculations of the special functions were performed with sufficient accuracy to ensure that T > S + L in all cases.) The value of E L limits the size of the largest of the (0, 0) terms, namely G 0,0 E , to a value on the order of 10 3 times the exact Green's function G (denoted ''G Exact'' in the tables), as expected. The numerically exact Green's function has been calculated using a pure spectral method, shown in (3).
[23] The second tables shown for each case (part (b)) show the values of G 0,0 E andG 0,0 E using E L obtained for        different values of L, keeping the frequency fixed fairly high such that d = a = 5.5 l 0 . It can be seen that the largest of the (0, 0) terms, G 0,0 E , has a magnitude on the order of 10 L times the magnitude of the total Green's function, as expected.
[24] If the number of significant digits desired for convergence S spec is specified, the calculated summation limits for the spectral series, P cal and Q cal , and the calculated summation limits for the spatial series, M cal and N cal , can be calculated using the formulae derived previously. These values are shown in part (c) of the tables for each case, along with the actual values P act , Q act and M act , N act needed to achieve convergence to S spec significant figures, obtained numerically. Also shown is S act , the actual number of significant digits to which the Ewald method has converged using the formula where G tot E = G E +G E is the value obtained from the Ewald method after summing the two series using P cal , Q cal , M cal and N cal . The frequency is again fixed such that d = a = 5.5 l 0 . For the spectral and the spatial series, the adjustment factor b = 1 works in all cases but is excessively conservative. A factor of b = 2 was assumed for the 1-D cases, and b = 4 for the 2-D cases. The agreement between the actual and specified values of S is good, especially for larger values of S, with S act > S spec in all cases except one (the 2-D case with S spec = 5).

Conclusions
[25] The Ewald method is a very efficient method for calculating the periodic free-space Green's function for three different cases: a 1-D array of line sources, a 1-D array of point sources, and a 2-D array of point sources. However, as noted in Kustepeli and Martin [2000] the method suffers from accuracy problems at high frequency due to a loss of significant figures that occurs from a cancellation when combining the spectral and spatial series that appear in the method. The method proposed here determines the ''best'' value of the splitting parameter E that appears in the method to yield the fastest convergence of the Ewald sums while limiting the number of lost significant digits to a specified level L. The derivation has been presented in a unified fashion so that all three cases are treated together, with Table 1 giving the parameters needed for the different cases. The predicted loss of significant digits is verified through numerical simulations and the results illustrate the accuracy of the proposed formula.
[26] Approximate expressions for the summation limits required to achieve a specified convergence accuracy to S significant figures for the Green's function were also formulated and tested for the three different cases. In these expressions the value of S is arbitrary, except that it should satisfy the constraint that S < T À L, where T is the number of digits in the arithmetic. Again, a unified derivation has been presented, with Table 2 summarizing the parameters to be used for each of the three different cases. The specified number of significant digits S desired for convergence was compared to the actual number of significant digits of the resulting series and found to be in good agreement in almost all cases, thereby validating the formulas.

Appendix A
[27] In this appendix we summarize the asymptotic evaluation of the exponential-integral series that appear in the 1-D point-source spectral series and the 1-D line source spatial series. The asymptotic evaluation of the terms involving the complementary error function that appear in the 1-D point-source spatial series, the 1-D line source spectral series, and the 2-D point-source spatial and spectral series, uses (10) and is straightforward.
[28] In performing the asymptotic evaluations of the exponential-integral series that appear in (5) and (6), the series have been cast back into their integral forms, namely the forms shown in equation (15) in Capolino et al. [2007] and equation (11) in Capolino et al. [2005]. These are the forms appearing in (11) and (12), respectively. In the following, case 1 denotes the 1-D line source array and case 2 denotes the 1-D point-source array. For each of these two cases, part (a) is defined as the subcase involving the determination of the splitting parameter, whereas part (b) is the subcase involving the determination of the summation limit. To determine the splitting parameter, either k or k r0 is the large parameter in the asymptotic evaluation. To determine the summation limits, either R m or the wave number k rp is the large parameter in the asymptotic evaluation. In (11), we made the substitution s = -u for case 2a. For case 2b we then used the notational change s = u for consistency. (Note that in (11) the upper limit denotes the complex point at infinity, which has been taken as positive infinity for the derivation in case 2b and negative infinity for the derivation in case 2a).
[29] In all four cases, each of the integral expressions can be put into the form where W is a real positive parameter that becomes large. The terms s 0 , W, and the functions f(s) and g(s) are defined in Table A1 for each of the four cases. For large W we use the integration-by-parts asymptotic approximation [Felsen and Marcuvitz, 1994;Bleistein and Handelsman, 1986] to obtain Table A1 summarizes for each of the four cases the integral to be evaluated (the first row), the final asymptotic results (the last row), and all of the parameters that appear in (A1).