Eﬃcient computation of the 3D Green’s function for the Helmholtz operator for a linear array of point sources using the Ewald method q

The Ewald method is applied to accelerate the evaluation of the Green’s function (GF) of an inﬁnite equispaced linear array of point sources with linear phasing. Only a few terms are needed to evaluate Ewald sums, which are cast in terms of error functions and exponential integrals, to high accuracy. It is shown analytically that the choice of the standard ‘‘opti-mal’’ Ewald splitting parameter E 0 causes overﬂow errors at high frequencies (period large compared to the wavelength), and convergence rates are analyzed. A recipe for selecting the Ewald splitting parameter is provided.


Introduction
In applying numerical full wave methods like the Method of Moments (MoM) or Boundary Integral Equations (BIE) to periodic structures made of conducting or dielectric electromagnetic scatterers, fast and accurate means for evaluating the periodic Green's function (GF) are often needed, in both microwave and optical frequencies where also other techniques could benefit by the fast evaluation of the periodic GF. The GF can also be used directly in analogous acoustic problems.
Among various techniques to accelerate computation of the periodic GF is the Ewald method, originally developed by Ewald in [1] to efficiently calculate the electrostatic scalar potential of a 3D periodic distribution of point charges. The method has been applied to a large variety of problems ranging from solid state physics and chemistry to engineering.
In the following we describe the application of the method to the electrodynamic or acoustic problems such as the solution of the Helmholtz equation with periodic boundary conditions. We list here the Ewald methods that have been previously developed, and identify our own contribution.
Jordan et al. extended the Ewald method to the case of the free space GF for three-dimensional (3D) problems with 2D periodicity (i.e., a planar array of point sources) in [2]. Its application in evaluating GFs for multilayered media is treated in [3][4][5][6]. The Ewald method is extended in [7][8][9] to 2D problems with 1D periodicity (i.e., a planar array of line sources). An analogous treatment of the GF for an array of line sources is reported in [10] where a more general view is also addressed, as well as in [11] where Schlö milch series that may arise in diffraction theory are efficiently evaluated. We also point out that in [8] a comparison is made between various representations of the periodic GF concluding that the Ewald method is one of the most advantageous. For a mathematical review of the Ewald method in general, see [12]. For the case of a planar 2D-periodic array of point sources, the critical distance from the array plane beyond which the Ewald method is not advantageous compared to the standard spectral grating lobe series is analyzed in [13].
The evaluation of the GF for a rectangular cavity using the Ewald method is reported in [14], and its dyadic form in [15]. An application of the Ewald representation to hybrid methods involving a boundary integral equation and finite element method, for periodic structures, is given in [16]. In [17] the Ewald method is combined with another very effective method to analyze scattering by periodically etched layered conducting surfaces such as are commonly used for frequency selective surfaces.
Here we accelerate for the first time with the Ewald method the GF pertaining to a linear array of point sources (i.e., a 3D problem with 1D periodicity). The formula representations obtained in the above mentioned studies, i.e., the GF for geometries such as the 2D-periodic array of point sources [2], and the 1D-periodic array of line sources [7][8][9][10], cannot be applied to this new geometry. The standard purely spectral sum representation for the GF cannot be used for this new particular configuration because it diverges when the observation point lies on the axis of the linear array, thus further proving the utility of the Ewald method, which exhibits a Gaussian convergence. A numerical analysis of the convergence rate is also provided, supported by analytical considerations. As previously noted in [18], at high frequencies such that the wavelength is somewhat smaller than the period, computing the Ewald series using the optimal splitting parameter E 0 (see [2]) may yield inaccurate results due to the finiteness of machine precision. In [18] it is suggested that such numerical inaccuracies may be improved by increasing the value of the splitting parameter E as frequency increases, but no guidelines for choosing the parameter are given. Here we present an algorithm for choosing the Ewald splitting parameter E that extends the efficiency of the method when the wavelength is somewhat smaller than the periodicity. The proposed algorithm is efficiently applied to periodic structures when the observation point is near the axis of the linear array of point sources.

Statement of the problem
For several numerical methods such as MoM or BIE applied to periodic structures with period d along z, it may be convenient to find the scalar potential radiated by an infinite linear array of point sources located at r 0 n ¼ n dẑ, where bold symbols define vector quantities, the caretˆdefines unit vectors, and n = 0, ± 1, ± 2,. . . The point sources may be excited with a linear phase progression e Àjk z0 z , where k z0 is the longitudinal wavenumber, that may originate by an impinging plane wave (impressed wavenumber) or by generators within each cell with a linear phase progression. An e jxt time dependence, with j ¼ ffiffiffiffiffiffi ffi À1 p , is assumed and suppressed throughout this paper. If one adopts an e Àixt time dependence, with i ¼ ffiffiffiffiffiffi ffi À1 p , all intermediate and final formula here derived apply as well after complex conjugation.
The GF for such a problem periodic along z, with period d, is given by the outgoing solution of the Helmholtz equation $ 2 G(r, r 0 ) + k 2 G(r,r 0 ) = À d p (r 0 ), where the periodic (along z) forcing term d p (r 0 ) is written as d p ðr 0 Þ ¼ P 1 n¼0 dðx 0 Þdðy 0 Þdðz 0 À ndÞe Àjk z0 nd . The GF solution of the Helmholtz equation may be represented as the spatial sum of spherical waves where R n ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q 2 þ ðz À ndÞ 2 q is the distance between the observation point r qq þ zẑ and the nth source point r 0 n ¼ n dẑ (see Fig. 1). Because of the symmetry of the problem, we use cylindrical coordinates and denote an arbitrary point as r = (q, z). Terms in the series (1) are of order n À1 exp(Àjnh), where h = (k z0 + k)d for large positive n, and h = (k z0 À k)d for large negative n, so that the series is extremely slowly convergent. The series diverges for h = 2mp, where m is an arbitrary integer, because the spherical contributions arising from at least one side of the linear array would all add in phase. In other words, one of the Floquet waves (see what follows) is propagating at grazing.
An alternative spectral series representation of this GF may be obtained by applying to (1) the infinite Poisson sum formula [19, p. 117], is its Fourier transform. This yields the GF representation as a sum of cylindrical harmonics, but it is also slowly convergent. Here H ð2Þ 0 is the Hankel function of zeroth order and second kind, is the qth Floquet wavenumber along z, and is the qth transverse Floquet wavenumber, along q. Convergence of the field far away from the linear array requires that Im k qq < 0, as can be inferred from the large argument expansion of the Hankel function. For a real excitation wavenumber k z0 and for radiation in a lossless medium (real k), q-harmonics such that k qq > 0, for k zq < k, index radially propagating Floquet waves, while when k qq ¼ Àj½k 2 zq À k 2 1=2 , for k zq > k, q indexes radially attenuating Floquet waves.
Note that the expression (4) has the important drawback that it cannot be used when the observation point lies on the line of the array, namely when q = 0, a situation that often happens in the MoM or BIE. Here we provide a representation of the GF (1) or (4) as a rapidly converging series, that can be evaluated also when q = 0.  where the integration path in the complex space is shown in Fig. 2. For simplicity the homogeneous medium is assumed to have small losses and the complex ambient wavenumber k = k r + jk i has k i < 0. The particular shape of the integration path is dictated by constraints for the integrability of the integrand of (7), which is explained as follows. It is useful to write the ambient wavenumber in polar coordinates k ¼ jkje j/ k , with / k = arctan(k i /k r ) < 0, and s = s r + js i = jsjexp(j/ s ) with / s = arctan(s i /s r ). Indeed, convergence of (7) is ensured by requiring that Re[k 2 /s 2 ] < 0 for s ! 0, and Re½R 2 n s 2 > 0 for s ! 1, i.e., 3 4 p

The Ewald Green's function transformation
which are satisfied by the path in Fig. 2. The Ewald method is obtained by splitting the integral in (7) into two parts, The spatial term G spatial exhibits a Gaussian exponential convergence in n as can be inferred by noticing that ReðsÞ > E on the complex integration path and for large n the integral in (12) behaves like 1 ffiffi is the complementary error function. For large argument it is approximated by erfcðzÞ $ e Àz 2 =ð ffiffiffi p p zÞ, which highlights the Gaussian n-convergence of (12). The spectral term (11) instead is slowly converging as 1/n, since the original series (1) converges as 1/n, and requires the use of the Poisson transformation to accelerate its convergence.

The spectral part of the Green's function
To accelerate the spectral part (11) of the GF in (10) we first apply the Poisson's sum formula, (cf., discussion before (2)), with Im s Re s 0 k k Fig. 2. Path of integration. Expressing k = jkjexp(j/ k ) (Imk 6 0, thus / k < 0), and s = jsjexp(j/ s ), the region of convergence of (7) is given in (8) and (9).
The evaluation of the Fourier transformf ð2pq=dÞ, defined in (3), is carried out by first changing the order of integration, which is possible since Re(s 2 ) > 0 on the integration path shown in Fig. 2. More details on these mathematical steps are in [9, Section IV-A] where an analogous problem is treated. Then, the change of variable (z 0 , s) ! (u, s), with u = (z À z 0 )s, after evaluation of the resulting u-integral, leads tõ Note that Reðk 2 qq Þ, with k 2 qq ¼ k 2 À k 2 zq , is negative for large q resulting in evanescent Floquet modes, hence the above integral is exponentially converging for large q. Let u ¼ 1 4s 2 , with ds s ¼ À 1 2 du u , so that the above integral is rewritten as where the mapped integration path goes to infinity in the region between p 2 and 3p 2 for real k qq (propagating modes), and in the region between À p 2 and p 2 for imaginary k qq (evanescent modes). The same argument applies when small ambient losses are present, while a deeper analysis is required for complex propagation constants k zq .
Next, following a procedure analogous to [7], [9, Section IV-B], we utilize the Taylor expansion e Àq 2 =ð4uÞ ¼ P 1 p¼0 ðÀ1Þ p ðq=2Þ 2p =ðp!u p Þ. After using the change of variable w ¼ 4E 2 u, an integration of the series term by term leads tõ in which E p (x) is the pth order exponential integral defined by [20, p. 228] and by the recurrence relation [20, p. 229] The path of integration in (18) should exclude the origin and not cross the negative real axis [20, p. 229]. For evanescent waves, Àk 2 qq > 0 and one can also use the definition E p ðxÞ ¼ R 1 1 e Àxt t Àp dt, valid for Re x > 0. Substituting (17) into the Poisson transform of the original expression (11) leads to G spectral ðr; r 0 Þ ¼ 1 4pd The exponential integral E 1 (x) may be evaluated using the algorithm of [20, Section 5.  (20) is rapidly convergent. Indeed, only a handful of q terms and around ten p terms are necessary to reach a good accuracy, as shown in Section 5.

Choice of the branch
The exponential integral (18) has a branch cut along the negative real axis. Note that in absence of ambient losses, and for a real phasing wavenumber k z0 along the array, all propagating Floquet waves have k 2 qq ¼ k 2 À k 2 zq > 0 and thus the argument of the exponential integral in (20) is negative real, resulting in an ambiguous branch condition.
To resolve the ambiguity we imagine small ambient losses such that ImðÀk 2 qq Þ > 0 and therefore the branch can be chosen automatically and one can use (18) directly. Though the assumption of losses is useful for the branch determination, in lossless media one may prefer to evaluate E p+1 (x) obtained using the formula E 1 (Àx + j0 + ) = À Ei(x) À jp, [20, p. 228], with where P denotes principal-value integration. For small argument the exponential integral Ei can be evaluated as where c % 0.57721566 is the Euler constant. Note that once E 1 (Àx + j0 + ) is obtained, E p+1 (Àx + j0 + ) follows from (19).

The spatial part of the Green's function
The spatial part of the GF in (12) already exhibits Gaussian convergence; here we reexpress it in terms of complementary error functions. Consider a generic n-term of the spatial series in (12) With the change of variable of u = R n s + jk/(2s), or s ¼ ðu þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 À j2kR n p Þ=ð2R n Þ, (22) is transformed to Now, the second integral is transformed by using the change of variable w ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 , and u 2 = w 2 + j2kR n , yielding Inserting this result into (12) leads finally to G spatial ðr; r 0 Þ ¼ 1 8p where the complementary error function erfc is defined in (13). For its fast evaluation, see [21].

The singular spatial contribution
The GF often appears convolved with induced or equivalent currents within a reference cell of the periodic structure (which we take to be the n = 0 cell), and special treatment is needed for the case when the observation point r is close to the source point in the cell, resulting in a vanishing R 0 . From the GF representation in (1) we know that when R 0 % 0, one has G(r, r 0 ) % 1/(4pR 0 ). This singular term is fully contained in the n = 0 term in the spatial sum in (25), since erfcðR n E AE j k 2 EÞ % erfcðAEj k 2 EÞ, and erfcðþj k 2 EÞ þ erfcðÀj k 2 EÞ ¼ 2 [20].

Asymptotic convergence of series in G spectral and G spatial
The convergence properties of G spectral and G spatial , in (20) and (25), respectively, are analyzed here. First we derive the so called ''optimum'' Ewald splitting parameter E 0 , and then we show that in some specific important cases this choice leads to a high-frequency breakdown.

The optimum Ewald splitting parameter E 0
In the spectral sum G spectral , for large q, we can approximate k zq % 2pjqj/d and k qq % Àj(2pjqj/d). By using the large argument asymptotic expansion for the exponential integral E p+1 (z) $ e Àz /z, the qth term of the series in (20) is asymptotically approximated for large q as $ e Àj2pqz=d e ÀðqEÞ 2 4pd which exhibits Gaussian q-convergence. In the spatial sum G spatial , for large n, we can approximate R n % jnjd, and the two error functions in G spatial are approximated by their asymptotic expansion for large argument as in the text after (13). Thus, the nth term of the series in (25) which exhibits Gaussian n-convergence. When truncating the two series in (20) and (25) as P Q q¼ÀQ and P N n¼ÀN , respectively, resulting in a total number of terms 2Q + 1 and 2N + 1, the optimum Ewald splitting parameter E 0 , minimizes the total number (2Q + 1) + (2N + 1) of necessary q, n-terms.
The total number of terms needed is determined by observing that to obtain roughly the same number of significant digits of accuracy in each sum in (20) and (25), the asymptotic approximation of their terms in (26) and (27), respectively, must have the same exponential factor e Àr 2 , and thus we set ðQpÞ 2 =ðEdÞ 2 ¼ ðN EdÞ 2 r 2 . The minimum number of terms needed to achieve this accuracy is thus M tot ¼ 2r½1=ðEdÞ þ Ed=p þ 2. The optimum E 0 parameter is obtained by imposing oM tot =oE ¼ 0, which leads to We note that the choice of the optimum splitting parameter E ¼ E 0 results in both series, G spectral in (20) and G spatial in (25), converging asymptotically with identical Gaussian convergence rates, $exp(Àpq 2 ) and $exp(Àpn 2 ), respectively.

The high-frequency breakdown of the Ewald representation
At high frequencies (large k), some of the leading terms in the Ewald representation in both series (20) and (25) can become very large but nearly cancelling, not only causing loss of accuracy but also leading to overflow errors. In particular, exponentially large terms arise from the low-q terms in the spectral sum G spectral in (20) that correspond to radially propagating Floquet waves (those with real k qq ). Without loss of generality we suppose that the array is phased with a real k z0 and that it is in a lossless free space. The q = 0 term has the largest value of k qq , as is often the case. The exponential integrals related to the q = 0 term can be approximated as which causes numerical instability when the argument k 2 q0 =ð4E 2 Þ is very large. These large numbers tend to cancel a large number arising from the spatial G spatial in (25), with a resulting loss of significant figures. The largest counterpart terms in the spatial sum G spatial in (25) is provided by the two erfc functions in (25) that asymptotically behave like The appropriate index n in (30) is that which minimizes R n (often n = 0) because the large growth e k 2 =ð4E 2 Þ is compensated by e ÀR 2 n E 2 when R n is large. It should be noticed that this numerical instability arises from the Ewald splitting and is not physical. The sum of the large numbers arising from the spectral and spatial splitting of course analytically leads to the correct results, but numerically may result in loss of accuracy or overflow. Thus at high frequencies, we have a loss of precision due to cancellation of large numbers when summing G spectral and G spatial in (10). The cancellation problem may be avoided by requiring that k 2 q0 =ð4E 2 Þ < H 2 in the spectral sum, and that k 2 =ð4E 2 Þ À R 2 0 E 2 < H 2 in the spatial sum, where H 2 is the maximum exponent permitted. This leads to Under the approximations k z0 ( k and kR 0 ( H 2 we have E spectral % E HF % E spatial . The fact that E spectral and E spatial have different expressions suggests that an algorithm where these large numbers are numerically stabilized a priori (mixing spectral and spatial terms) is not straightforward, and is considered beyond the scope of this work. If one desires to use the optimum parameter E 0 in (28) (that minimizes the number of terms) with the constraint that the largest computed values in G spectral and G spatial are smaller than e H 2 that may cause loss of accuracy, one should also impose that E 0 P E HF . For example, if we require H 2 = 6 we need E > E HF % k=4:9 % 1:28=k. In other words, for H 2 = 6, the choice of the optimum parameter E ¼ E 0 is a good choice if d < 1.38k.
Even though the exponential growth is limited by the choice E > E HF , in order to have a computationally efficient algorithm we also require the p-sum in (20) to be rapidly convergent.  (20) relies on having ðqEÞ 2p =ðp þ 1Þ! negligible for p P P, i.e., ðqEÞ 2P =ðP þ 1Þ! < , with the desired error and P the number of p-terms necessary to achieve convergence (typically [10][11][12][13][14][15]. This implies that For example, if we require = 10 À6 , and P = 15, we have E < E Taylor % 1:75=q that has to be considered jointly with (28). In other words, the optimum parameter in (28) can be used if E 0 < E Taylor which is, for the numerical example just given, equivalent to imposing that the Ewald sum be used up to q < 1:75d= ffiffiffi p p % d. Thus for high frequencies, or equivalently for large interelement spacings d > k, the constraint (32) forces a choice of E other than ''optimum'' (28). We therefore suggest choosing for the E parameter in (20) and (25) whose expressions can be used for the GF evaluation up to a distance q determined by (33). When the frequency is sufficiently high that the Ewald parameter has to be chosen as E ¼ E HF , one should use the Ewald algorithm under the condition E HF < E Taylor , which is equivalent to q/k < [(P + 1)!] 1/2P H/p that may further restrict the radial distance q where the Ewald sum is accurate.

Numerical results: convergence
Various examples describe the convergence behavior of the two sums (20) and (25), for G spectral and G spatial , respectively, of the Ewald splitting (10). In all cases that follow, the linear array is phased, i.e., k z0 = 0.1k and the n = 0 point source is at (q 0 , z 0 ) = (0, 0).
In Fig. 3 we analyze the convergence rate of the Ewald sums using the percent relative error defined as where G exact is the GF reference solution evaluated via its spectral counterpart (4) with a sufficiently large number of terms to achieve accuracy up to eight decimal digits, and G Ewald is the same GF evaluated using the Ewald splitting (10) with (20) and (25). For the case q = 0, the spectral series (4) cannot be used, and the reference solution is evaluated with 10 8 terms of the spatial series (1). The percentage relative error is plotted versus summation limit parameters ±Q and ±N employed in sums (20) and (25), respectively (in particular, ÀQ 6 q 6 Q, and ÀN 6 n 6 N), resulting in a total number of terms of 2P + 1 and 2N + 1. The three curves are related to three observation points at q = 0, 0.01d, 0.1d, all with z = 0.1d, where d is the period, which is set to d = 0.05 k and d = 0.5k, with k = 2p/k the free space wavelength. The Ewald splitting parameter is chosen as E ¼ E 0 given in (28). A large total number P = 20 of p-terms in (20) has been used because we emphasize here convergence issues related only to N and Q. Fig. 3 shows that just a few terms, Q = 1 and N = 1, provide six digits of accuracy. The relative error cannot be further decreased by augmenting the number of terms Q, N in G spectral , G spatial because of the accuracy limits of the numerical subroutines that evaluate the error functions in G spatial and the exponential integral E 1 (x) in G spectral (see [21] for erfc(x), and [20, Sections 5.1.53 and 5.1.56] for E 1 (x)). Note that in both cases in Fig. 3 the period d is smaller than the wavelength, and that accepting exponential factors of the order of e H 2 with H 2 = 9.6, i.e., with roughly six digits of accuracy, implies E HF ¼ k=6:2 ¼ 1:01=k. According to (28), for d = 0.05k and d = 0.5k, the optimal Ewald parameter is E 0 ¼ 35:4=k and E 0 ¼ 3:54=k, respectively, and thus in both cases we have E 0 > E HF justifying the use of the optimum Ewald parameter E 0 .
To gain an idea of the exceptional convergence rate of the Ewald representation, independent of the accuracy error of the evaluation of the special functions E 1 (x) and erfc(x), Fig. 4 shows the rate of convergence of the two individual series G spectral and G spatial in (20) and (25), respectively, evaluated at the location r = (q, z) = (0.1,0.1)d (we have checked that the case with z = 0.1d would lead to analogous results). The relative error is defined as Err spectral;spatial ¼ jG Ewald;exact spectral;spatial À G Ewald;P ;M spectral;spatial j jG Ewald;exact spectral;spatial j Â 100 ð36Þ where G Ewald;exact spectral;spatial is either G spectral or G spatial evaluated with a sufficient number of terms to achieve high numerical accuracy. The error is plotted versus the summation limit parameters ±Q and ±N employed in the sums in G spectral and G spatial , respectively. The array is the same as that in Fig. 3(b), i.e., d = 0.5k. The convergence rate is shown for three different choices of the Ewald splitting parameter E. Note that the optimum parameter E ¼ E 0 leads to the same convergence rate for both G spectral and in G spatial (Fig. 4(b)), and the two series need the same number of terms to achieve the same accuracy. It is remarkable that so few terms are required to achieve high accuracy. For E < E 0 the spectral sum in G spectral converges more rapidly than the sum in G spatial . Vice versa for E > E 0 , as predicted by the asymptotic expansion of the spectral and spatial sums in (26) and (27), respectively.
To determine the number P of necessary p-terms used in the Taylor expansion in (17), the percentage relative error of the spectral part (20) of the Ewald representation is plotted in Fig. 5(a) versus the radial distance q of the observation point from the array axis. The number of Ewald terms in (20) is kept equal to Q = 2, which is enough to guarantee good accuracy as shown in Fig. 3. The relative error in Fig. 5(a) is defined using as reference solution the same expression (20) evaluated with a large number P = 60 of p-terms. A larger number P of terms is necessary for larger distances q from the array axis. We observe that in Fig. 5(a), P = 20 is satisfactory for good accuracy such as 10 À4 % until q = 1.3d. In Fig. 5(b), for the same settings of Fig. 5(a), we plot instead the percentage relative error of (10) with (20) and (25) (with Q = N = 2), taking as reference the spectral representation (4) of the GF, which is not affected by the accuracy limits of the numerical subroutines that evaluate the exponential integral E 1 (x) in (20) and the error functions erfc(x) in (25). This is why in Fig. 5(b) the accuracy never falls below a certain value (determined by the accuracy limits of erfc(x) and E 1 (x)) even for small q. Note that for P = 20 the accuracy becomes 10 À4 % at q % 1.3d. The trend is in accordance with the reasoning that led to (33). Repeating that discussion with P = 20, with a truncation error = 10 À6 , we would have E Taylor ¼ 2:2=q, and the condition (33) is equivalent to q < 1.24d.

Extension to high frequency
In Fig. 6(a) the frequency is further increased such that d = 5.5k and the high frequency breakdown problem described in Section 4.2 occurs. In this case we cannot choose the optimum parameter E ¼ E 0 . Indeed, if we did, we would have k=ð2E 0 Þ ¼ ffiffiffi p p d=k ¼ 9:75 in (25), yielding terms on the order of e 95 that would cause overflow or at least loss of accuracy. The Ewald splitting parameter E must therefore be chosen according to (34). The rate of convergence of the two series in (20), (25) is then different, and more terms (Q > N) in G spectral (20) are needed in order to maintain the same relative error in both (20) and (25). The convergence rates of the spectral and spatial sums of the Ewald representation are reported in (26) and (27), respectively. In Fig. 6(a) the error is plotted versus N, with Q = 4N. Note that again a few terms are necessary (N = 2 and Q = 8) to provide a good accuracy. The error cannot further decrease than that for N = 3 and Q = 12 by augmenting the number of terms Q, N in G spectral , G spatial because of the accuracy limits of the numerical subroutines that evaluate the error functions and the exponential integral. Fig. 6(b) shows the rate of convergence of the two individual series G spectral and G spatial in (20) and (25), respectively, evaluated at the location (q, z) = (0.1,0.1)d. The relative error is defined as in (36). The two series in G spectral and G spatial need different numbers of terms to achieve the same accuracy since we have chosen E ¼ E HF > E 0 . The spatial series G spatial needs only the n = 0 term and it is accurate enough to be off the scale of the plot; hence only the number Q of terms in the spectral series G spectral is given. The error in Fig. 6(b) does not decrease significantly before Q = 6 because for this case (d = 5.5k) spectral terms with jqj 6 Q = 6 are propagating and must be accounted for to obtain an accurate evaluation of the spectral part of the Ewald representation.
Note that the larger the Ewald parameter, the greater the number of terms needed to achieve the same accuracy. The three cases with E ¼ E HF , are relative to E HF ¼ 1:05=k ¼ 3:28E 0 , E HF ¼ 1:19=k ¼ 3:72E 0 and E HF ¼ 1:57=k ¼ 4:9E 0 , where the optimum parameter (not so optimum in this case) is E 0 ¼ 0:32=k. According to (34), these cases correspond to choosing H = 3, 2.64 and 2, respectively. Clearly, as described in Section 4.2, choosing a larger H (i.e., smaller E HF =E 0 ) bounds the growth of the number of terms in G spectral . Unfortunately, H cannot be chosen too large because, as explained in Section 4.2, it would generate values on the order of e H 2 . It should also be noted that choosing E ¼ E HF , with decreasing H = 3, 2.64, 2, automatically implies that the maximum values in the series G spectral and G spatial become smaller. For an array with d = 5.5k, values of H > 3 (i.e., E HF =E 0 < 3:28) would already cause a small loss of accuracy, while H < 2 (i.e., E HF =E 0 > 4:9) would require too many terms in the G spectral sum.

Conclusion
The Green's function (GF) for a linear array of point sources linearly phased has been expressed in terms of two series that exhibit Gaussian convergence as shown in Section 4.1. A description of the convergence properties, the analysis of the ''high frequency breakdown'' of the Ewald representation, and its remedy, have been analyzed theoretically and with the help of numerical experiments.
It is important to note that the purely spectral representation of the GF (4), besides being much slower than the Ewald representation, cannot be evaluated for q = 0, which is a very typical situation encountered when analyzing periodic structures of 3D elements with 1D periodicity. In this particular case (q = 0) the spatial representation in (1) (sum of spherical waves, which is in general even slower than the spectral representation in (4)) 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 well as the spectral method for q = 0), while the Ewald method can be evaluated also in this case. The Ewald representation is particularly useful for small q, while for large q (of the order of a period), one can use the well known spectral GF in (4) that is rapidly converging in this particular situation. Automatic optimal switching between these two representation for a given accuracy is beyond the scope of this work and will be investigated in the future, as well as a more thorough automatic algorithm for deciding the number of terms to be used in the series and the switching parameters, as has been done in [22] for the Ewald representation discovered in [2].  (20) and (25), respectively. One term (n = 0) in G spatial is sufficient for excellent accuracy. In this case the choice of the optimum parameter E ¼ E 0 would cause numerical errors, as explained in Section 4.2. Three different choices of the E parameters are analyzed here. Increasing the ratio E=E 0 results in a larger number of terms in G spectral to reach convergence. A smaller ratio E=E 0 would decrease the number of necessary terms in G spectral but it would result in a loss of accuracy.