Efficient computation of the 3D Green's function with one dimensional periodicity using the Ewald method

The Ewald method is applied to accelerate the evaluation of the Green's function of an infinite periodic phased linear array of point sources. Only a few terms are needed to evaluate Ewald sums, which are cast in terms of error functions and exponential integrals, to high accuracy


Introduction
In applying numerical full wave methods like the Method of Moments (MoM) to periodic structures, fast and accurate means for evaluating the periodic Green's function are often needed. Among various techniques to accelerate computation of the Green's function is the Ewald method, originally developed by P. P. Ewald in [1], and extended to the case of the free space Green's function for three dimensional (3D) problems with 2D periodicity (i.e., a planar array of dipoles) in [2]. Its application in evaluating Green's functions for multilayered media is treated in [3] and [4]. The Ewald method is extended in [5], [6], [7] to 2D problems with 1D periodicity (i.e., a planar array of line sources), while its application in evaluating Green's functions for a rectangular cavity is reported in [8]. Here we accelerate for the first time with the Ewald method the Green's function pertaining to a linear array of dipoles.

Statement of the Problem
The scalar potential radiated by the infinite linear array, linearly phased along z, is G (r, r ) = ∞ n=−∞ e −jkz0nd e −jkRn /(4πR n ) where k z0 = kd cos θ is the phase gradient along the array, d is the lattice spacing, θ is the radiation pointing angle, k is the free space wavenumber, and R n = ρ 2 + (z − z − nd) 2 , with ρ = x 2 + y 2 , is the distance between the observation point r ≡ (ρ, z) and the nth source point r n ≡ (ρ = 0, z = nd) (see Fig. 1). An e jωt time dependence is assumed and suppressed throughout this summary. Terms in the spatial series representation of G are of order 1/n for large n, so that the series is extremely slowly convergent. An alternative spectral series representation of this Green's function, in terms of cylindrical waves, also exists, but it is also slowly convergent. Here H zq is the transverse Floquet wavenumber, along ρ. Convergence of the field far away from the linear array requires that Imk ρq < 0, as can be inferred by using the large argument expansion of the Hankel function. Note that the spectral representation (1) has the important drawback that it cannot be used when the observation point lies on the line of the array, namely when ρ = 0, a situation that often happens in the MoM.
Here we provide a representation of the Green's function as a rapidly converging series, that can be evaluated also when ρ = 0.

The Ewald Green's Function Transformation
Following [2] the Green's function for a single source is represented as The Ewald method method consists of splitting the integral in (2) into two parts, which also determines the splitting G r, r = G spectral r, r + G spatial r, r The Ewald splitting parameter E has to have a proper value to make the splitting efficient. The spatial term is exponentially converging in n while the spectral term is slowly converging and requires the use of the Poisson transformation to accelerate convergence. Skipping here the details of the derivations, similar to hat presented in [7], closed form formulas are obtained for the spectral and spatial terms.
in which erfc is the complementary error function and E q (x) is the qth order exponential integral defined as [9, p. 228 This function has a branch cut along the negative real axis and thus defines the solution for propagating Floquet modes with a small loss and propagating modes for lossy medium. From a numerical point of view, only the exponential integral E 1 (x) needs to be evaluated numerically using, for instance, the algorithm of [9,Sec. 5.1]. Higher order exponential integrals may be evaluated by the recurrence relation E q+1 (x) = 1 q [e −x − xE q (x)], q = 1, 2, 3, .... Note that in absence of losses all propagating Floquet waves have k ρp = k 2 −k 2 zp > 0 and thus the argument of the exponential integral in (7) is negative, resulting in an ambiguous branch condition. To resolve the ambiguity, we imagine small ambient losses so that Im(−k 2 ρp ) > 0 and therefore the branch can be chosen automatically.
Numerical Results: Convergence We analyze in Fig. 2 the convergence rate of the Ewald sums using the percent relative error defined as where G exact is the Green's function reference solution evaluated via its spectral counterpart with sufficiently large number of terms to achieve accuracy up to seven decimal digits (2000 terms are sufficient), and G Ewald is the same Green's function evaluated using the Ewald splitting (3) with (6) and (7). The percentage relative error is plotted versus summation limit parameters ±N , ±P in sums (6) and (7) resulting in a total number of terms of 2N +1 and 2P +1, respectively. The array is phased, i.e., k z0 = 0.1k. The n = 0 point source is at (x , y , z ) = (0, 0, 0), and the two curves are related to two observation points at (x, y, z) = (0.01, 0, 0.1)λ 0 and (x, z) = (0.1, 0, 0.1)λ 0 , respectively, where λ 0 = 2π/k is the free space wavelength. The period is set to d = 0.6λ 0 , typical of array antennas. In Fig. 2 the error is plotted versus N , with P = N . The Ewald splitting parameter is chosen as E = √ π/d (see [7]). A large total number Q = 30 of q-terms in (7) has been used because we emphasize here convergence issues related only to N and P . Note the sum limits in Eqs. (6) and (7) are plus and minus N,P, and positive Q. In our cases the relative error cannot be further decreased by augmenting the number of terms P, N in G spectral , G spatial because of accuracy limits of the numerical subroutines that evaluate the error functions in G spatial and the exponential integral E 1 (z) in G spectral (see [10] for erfc(z), and [9, Secs. 5.1.53, 5.1.56] for E 1 (z)).

Conclusion
The Green's function for a linear array of point sources linearly phased has been expressed in terms of two series that exhibit gaussian convergence that can be demonstrated by following the steps in [7, Sec.V]. It is important to note that the spectral representation of the GF (1), besides being much slower, cannot be evaluated for ρ = 0, which is a very typical situation. In this particular case the even slower spatial representation (sum of spherical waves) can still be used, rendering the Ewald method even more desirable. When it is required to evaluate the Green's function for complex phasings, the spatial representation cannot be evaluated because it would diverge (as  (3) with (6) and (7) evaluated at two observation points at r = (x, y, z) = (0.01, 0, 0.1)λ0 and at r = (x, y, z) = (0.1, 0, 0.1)λ0. Percentage relative error versus number of terms N (P = N ) in the sums. The period is d = 0.6λ0, the n = 0 source is at r = (x , y , z ) = (0, 0, 0). well as the spectral method for ρ = 0), while the Ewald method can be evaluated also in this case.