Acoustic scattering from a one-dimensional array; Tail-end asymptotics for efficient evaluation of the quasi-periodic Green’s function

• Plane wave acoustic scattering from an infinite periodic array of cylindrical scatterers is analysed using a boundary integral equation approach. • A new asymptotic correction term for the quasi-periodic Green’s function is introduced. • Results of the new approach are compared against existing methods. • Transmission and reflection coefficients for a variety of scatterer cross-sections are calculated using the boundary element method, showing good agreement with previous results. Motivated by the problem of acoustic plane wave scattering from an infinite periodic array of cylindrical scatterers, we present a new and easily-implemented way of calculating the quasi-periodic Green’s function. This approach is based on an asymptotic expansion of the summand in the quasi-periodic Green’s function in order to derive a tail-end correction term, allowing for a rapid and accurate approximation of the function. The tail-end approximation is shown to have much better and faster convergence properties than the usual truncation approach and competes very well with state-of- the-art alternative techniques. This method is then combined with a boundary element scheme to calculate the transmission and reflection coefficients associated with arrays of cylinders of different cross-sections and varying aspect ratios. The results are validated against the existing literature and by independent finite element calculations. a more periodic fashion than the truncated sum alone, which varies in amplitude. The comparison of the errors shows a great deal of variation in the truncated sum alone, whereas the sum with correction becomes very small method of calculating the two-dimensional free-space quasi-periodic Green’s function. We began by formulating the boundary integral equation for the scattering of an acoustic plane wave from a one-dimensional array of identical cylinders of arbitrary cross-section, and derived the form of the quasi-periodic Green’s function. After truncating the infinite sum present in the quasi-periodic Green’s function, we then used an asymptotic expansion for the remaining tail ends to produce the correction terms. This approach provides a fast and efficient way to approximate the function in a manner which is straightforward to implement. We compared our new method of calculation for the quasi-periodic Green’s function against state-of-the-art methods, including Ewald’s method and Kummer’s transformation. We found that our results were in good agreement with the results of [11], and although


Introduction
The interaction of waves with a periodic array of scatterers is a canonical problem that can be traced back to Rayleigh and his work on the properties of gratings [1]. He was influenced by the work of Fresnel on the diffraction of light, but scattering by infinite arrays is a problem that has a very wide range of applications in various fields. For example, [2] considered the scattering of electromagnetic waves from a grating of cylinders, which has influenced the design of antennae. The scattering of waves from the periodic structures present in composite materials used in underwater acoustics is key to characterising their properties as seen in [3]. More recently, scattering from periodic arrays has been applied to the design of phononic crystals, an overview of which can be found in [4], whilst the diffraction of waves by an array of axisymmetric ice floes was used in [5] as a model for marginal ice zones. The rapidly expanding field of metamaterials includes many designs that utilise periodic arrays of resonators or holes; a good collection of references can be found in [6].
Given such a broad range of applications, it is unsurprising that numerous methods of solution to the problem of scattering from arrays are to be found in the literature. The separation of variables and multipole method for an array of N cylinders was originally introduced in [7] and later extended to the infinite case in [8]. Twersky also proposed a method of solution for a one-dimensional array of parallel cylinders using integral equation methods [9], though his paper lacks any implementation. Achenbach et al. were able to solve the integral equation formulation for an oblique incident wave on a one-dimensional array of cylinders by the use of the boundary element method (BEM) [10]. The integral equation method is a popular choice for the solution of exterior scattering problems; unlike finite element methods, it does not require the entire domain to be discretised or the use of perfectly matched layers, for example, to deal with the radiation condition at infinity. The method may also be applied to scatterers of arbitrary cross-section, unlike those methods that rely on a circular or elliptical cross section. The chief drawback of this approach, and the motivation for our study, is that in formulating the boundary integral equation, we are required to deal with a quasi-periodic Green's function -in the case of acoustic waves, a sum of Hankel functions, which is notoriously difficult to compute. In order for any BEM scheme to succeed, an efficient means of computing the quasi-periodic Green's function must be found.
Numerous approaches for improving the speed and accuracy of such computations can be found in the literature. Linton provided a review of the alternative methods for calculating the quasi-periodic Green's function for the Helmholtz equation in [11], including Ewald's method as proposed in [12], Kummer's transformation, integral equation methods as discussed in [13], and lattice sums, as used in [14] and [15]. In this paper, we present a new method to achieve efficient computation of the quasi-periodic Green's function. The scheme is easy to implement in the form of (3.4) using expressions given in (3.5) and shows good performance over a range of parameters, save certain 'resonant' frequencies for a particular incident angle, which we discuss in brief below. We anticipate that the method presented in this paper may also be extended to the fully three-dimensional problem of an array of axisymmetric scatterers.
We begin in Section 2 by setting out the boundary integral equation formulation of acoustic plane wave scattering from a one-dimensional array of infinite cylinders of arbitrary cross-section and find the form of the quasi-periodic Green's function. In the following section we seek to derive an efficient computational form of the quasi-periodic Green's function, and its normal derivative, by truncating the infinite sum and expanding the remaining 'tail ends' asymptotically, which results in a series of correction terms that can be used to compute the function with increasing accuracy. Numerical results demonstrate the efficacy of the correction terms in comparison to the truncated sum only. We then implement our newly derived method in a BEM scheme and calculate the reflection and transmission coefficients for arrays of cylinders with a variety of cross sections, which are then validated by comparison with FEM results. Furthermore, the results for the circular cylinders show good agreement with those of [10] and [16].
As should be expected from previous studies, the numerical results show poorer convergence in the region close to the cut-on frequency of a newly propagating mode. This phenomenon is known as Wood's anomaly, and corresponds to the frequency and angle of incidence at which one of the scattered wave modes propagates along the array, resulting in a resonance effect. Whilst an important issue, we choose not to deal with these anomalous frequencies here, and direct the reader to [17], where a method to deal with these resonant frequencies is developed. Unfortunately, these methods rely on exploiting the circular cross-section of the scatterers, and hence are not applicable to non-circular cylinders. An alternative method is also proposed in [18].

Problem formulation
We consider a periodic array of cylinders with arbitrary cross-section and axes aligned parallel to the x 3 axis, as shown in Fig. 1. We denote the mth cylinder by V m , with surface ∂V m . We seek the solution φ(x), with x = (x 1 , x 2 , x 3 ), to the scattering problem, such that is the time-harmonic acoustic pressure, and the incident field is an incoming time-harmonic plane wave, which has the form Here we have introduced the wavevector k = −k(sin ψ 0 cos θ 0 , sin ψ 0 sin θ 0 , cos ψ 0 ), defined in terms of the incident angles (ψ 0 , θ 0 ), as illustrated in Fig. 1. The total field satisfies the Helmholtz equation in the domain D exterior to the cylinders, and we impose either Dirichlet or Neumann boundary conditions on the surface of every cylinder, so that for every V m , where n is the unit normal to the boundary of the cross-sectional surface ∂V m pointing into the cylinder. The total acoustic field is decomposed into the sum of the incident field, φ in , and scattered field φ sc Given the translational invariance of the problem in the x 3 direction, we may rewrite the scattered field as It may then be shown that φ sc (x 1 , x 2 ) also satisfies the two-dimensional Helmholtz equation with modified wavenumber Having factored out all x 3 dependence, we may now consider the decomposition Hence, we need only solve the two-dimensional Helmholtz equation with reduced wavenumber k. From now on, for convenience, we write φ(x 1 , x 2 ) = φ(x 1 , x 2 ) and x = (x 1 , x 2 ). Assuming x ∈ D, exterior to the cylinders, it is possible to rewrite the problem as a boundary integral equation over the surfaces ∂V m of the cylinders. For the Dirichlet problem, this gives and for the Neumann problem The free-space Green's function G is here given by We can now exploit the quasi-periodicity of the domain to reduce (2.2) to an equation over a single reference cylinder, V 0 . For an array of cylinders with spacing d, we define the zeroth cell containing V 0 as for field point x 0 , and source point ξ 0 ∈ S 0 , the boundary of the cross-section of V 0 . We can then define periodic coordinates that map x 0 and ξ 0 to x p in the pth and ξ q in the qth cell respectively, i.e.
Using the periodic coordinates, we can then rewrite the incident field in terms of which suggests that the total field also has the form We may now rewrite the integral equation (2.1) in terms of the zeroth cell coordinates, starting by expressing After multiplying through by a factor of e ikpd sin ψ 0 cos θ 0 and setting m = p − q, we are then able to write (2.1) in terms of an integral equation around the cross-sectional boundary S 0 of the reference cylinder V 0 , i.e.
Similarly, for the Neumann case, where we are now using the quasi-periodic Green's function, defined as Whilst this approach is amenable to solution by BEM, it also creates problems numerically since the sum in (2.6) is slow to converge in the context of both the Dirichlet and Neumann problems. A number of methods have been proposed to speed up this convergence e.g. [11]. We now present a new, efficient way of evaluating this infinite sum quickly and accurately by determining asymptotic correction terms to the tail ends that remain after truncation of the sum in (2.6). This approach is extremely straight forward to implement once the corrections terms have been derived.

Asymptotic corrections for the quasi-periodic Green's function
In order to speed up computation of the quasi-periodic Green's function (2.6), we seek to truncate the sum at some integer M, and then approximate the remaining tail ends of the sum, i.e. the sum for |m| ≥ M, using an asymptotic expansion for the summand for each term |m| ≥ M. Let us define the truncated sum as and seek asymptotic corrections σ ± (M) such that where σ + and σ − denote the correction terms for the tail ends m ≥ M and −m ≤ −M respectively. In order to study the asymptotic behaviour of the summand, we will introduce variables u = x 0 The periodic distance (2.7) may then be written as r m = √ (u + dm) 2 + v 2 . We make the assumption that kd, ku and kv We first seek the asymptotic form of the Hankel function, H We can also expand (2.7) for large m to give Therefore, expanding e iz √ z for large z, and substituting z = kr m with m ≫ 1, we arrive at the following expansion for where we have defined with ± corresponding to positive and negative m respectively. The correction terms for the tail ends are therefore given It is possible to rewrite σ ± in a more computationally efficient manner by removing the infinite sums in (3.2), as we shall now outline. We first define α ± = kd (1 ± cos θ 0 ) , so that due to the form of the summand, the right hand side of (3.2) can be rewritten as a sum of terms of the form In this form, we note that the sums may be written in terms of the Lerch transcendent [20] which is defined as We may then exploit the fact that the Lerch transcendent has the integral representation [21] Λ[z, in combination with Watson's Lemma to find the asymptotic form of these integrals for large M, provided that z ̸ = 1, i.e. α ± is not an integer multiple of 2π (we will discuss the case z = 1 in more detail below). So, as M → ∞, ) .

(3.3)
We therefore rewrite σ ± in terms of Lerch transcendents, beginning with (3.2), ( Finally, recalling (3.1), we find the large M asymptotic form for G P is given by We can apply the same approach to the normal derivative of G P , for which the correction terms are defined as , and The above analysis will fail when we encounter values of kd for which α ± is an integer multiple of 2π . These values correspond to the cut-on frequencies of the array, when the system is known to exhibit either single or double resonance, as one or two modes begin to propagate along the array itself. Evaluation of the Green's function in these special instances has been discussed in [17], though only for circular cylinders.

Error analysis
We now wish to compare the performance of our approximation against existing methods and the accuracy of our results as we vary the truncation point M. [11] surveys numerous methods for computing the quasi-periodic Green's function, and compares their performance by computing the Green's function with each method up to a given tolerance, i.e. 10 significant figures. We compare his results with our approximation in Tables 1 and 2. The value of M given represents the minimum truncation point required to reach a stable answer of 10 significant figures in both the real and imaginary parts of the Green's function, as sought by Linton. It is clear that the truncated sum alone shows very poor convergence and little agreement with the other methods (5000 represents the maximum truncation point computed by [11]). The asymptotic correction shows good agreement with both Kummer's transformation and Ewald's method, but requires a larger number of terms to be computed in order to achieve this 10-figure accuracy. We note, however, that our method benefits from having only a single truncation parameter to modify, as opposed the to the four parameters that must be chosen for Ewald's method.
In order to compute the absolute and relative errors of our approximation, we first must have a reference value against which to compare. We use for this reference value G ref , which is found by taking a suitably large number of terms in the sum to ensure the value has converged to an accuracy of 10 −14 . We define the absolute error for the approximation with truncation point M as  Table 2 Table of values taken from [11] showing minimum M, or collection of parameters, required to achieve We can see in Fig. 2(a) that for the 4 values of k investigated, the absolute error does decrease like M − 7 2 , as we would expect from the form of (3.4).
Similarly, the relative errors shown in Fig. 2(b) also conform to the same power law. The plots also demonstrate the efficacy of the correction terms, with absolute and relative errors reaching O(10) −9 at a truncation point M = 100. These results are only valid away from the cut-on frequencies, at which point the Green's function diverges; close to these anomalous frequencies, many more terms are required to achieve a similar degree of accuracy. We can see in Fig. 3 that in order to reach 10 significant figure accuracy, a truncation point of M = 10 3 is sufficient for the majority of wavenumbers in the range k = (0, 2.5), except in the regions close to the cut-on frequencies, where a far greater number of terms is required to achieve the same level of accuracy. Other methods also struggle to compute accurately the quasi-periodic Green's function close to these resonant frequencies.
In most practical applications, this level of accuracy is not required; accuracy up to 5 significant figures is often considered to be sufficient. If we therefore relax the condition on the accuracy of our results from 10 to 5 significant figures, we can see from Fig. 4 that the required value of M is much lower, even close to the resonant frequencies. In such contexts that require a lesser degree of accuracy, our method represents an easily-implemented and efficient means to calculate the quasi-periodic Green's function.
Recalling the periodic behaviour of the Green's function in the x 1 direction, we can also compare to what extent our approximations satisfy this condition by taking the largest difference between any two peaks in the function over a range , say, as a measure of the error. In Fig. 5 we see that at M = 150, the value of the sum with the correction behaves in a more periodic fashion than the truncated sum alone, which varies in amplitude. The comparison of the errors shows a great deal of variation in the truncated sum alone, whereas the sum with correction becomes very small

Transmission and reflection results
We can now implement the correction term (3.4) in a BEM scheme, as outlined in the preceding section, in order to calculate reflection and transmission coefficients for infinite arrays of cylinders of different cross-sections. We first illustrate some results for the reflection and transmission coefficients, R n and T n , as defined in [10] for the Neumann boundary condition. These satisfy the power balance where the sum is taken only over the values of n for which the modes are non-evanescent. The coefficients α n are given in terms of the incident angle and spacing d by with α n = 0 giving the cut-on point of the nth mode. Given the appearance of α n in the denominator in (4.1), it is clear that the sum ceases to converge at these values. As can be seen in Fig. 6, close to the cut-on frequencies, we see spurious values above 1 introduced, which are a result of the numerical scheme breaking down -it is at these points that Wood anomalies appear, that is, the scattered wave propagates along the array creating a resonant effect. We provide these results for purposes of comparison with [10] and [16], showing good agreement away from the cut-on frequencies.
We note that the presence of the 'normalising factor' α 0 αn in the power balance used by Achenbach suggests that R n , T n do not in fact represent the physical quantities of reflected and transmitted energy in each mode. We therefore rescale our reflection and transmission coefficients by the 'normalising' factor to obtain so that R n , T n satisfy the power balance  In order to validate our results, we compare the values of the (normalised) reflection and transmission coefficients obtained from the BEM scheme with those calculated using the FEM software COMSOL Multiphysics, version 5.3. A sample of these results is shown in Fig. 7, where a truncation point of M = 600 was used. Quadratic isoparametric elements were used to discretise the boundary for the BEM scheme, and the number of elements, N = 48, was chosen to ensure numerical convergence within 10 −6 . Gaussian quadrature was used to numerically evaluate the integrals. As can be seen both the circular and elliptical cases show good agreement, though as expected both methods struggle close to the cut-on point. We can now compare the reflection and transmission for a variety of cross sections, as seen in Figs. 8 and 9. In Fig. 8, we compare cylinders with a 'pear'-shaped cross section to elliptical cylinders of a similar aspect ratio, to examine what difference the deformed shape has on the reflection and transmission of energy. It can be seen from Figs. 8(a) and 8(b) that the two cross sections give similar results, though there is some difference for the higher wave modes.
Similarly, we compare an array of upright ellipses in Fig. 9(a) to an array of cylinders with 'peanut'-shaped cross sections in Fig. 9(b) with the same aspect ratio. A far bigger difference can be seen, with the peanut-shaped scatterers showing much higher transmission and lower reflection than the elliptical cylinders. We expect that this is owing to the concave sides of the peanut shape allowing more of the incident wave to be transmitted through the array.
It is also interesting to compare the effects of aspect ratio and spacing of, for example, elliptical scatterers. In Fig. 10(a) we compare the transmitted energy of elliptical scatterers with different aspect ratios as the ratio between the major radius a and the period d increases. It is interesting to note the non-monotonicity of the curve for aspect ratio b a = 10, with a local maximum being reached close to a d = 0.4. Similarly in Fig. 10(b) we take scatterers with different normalised areas, and compare the transmission as the aspect ratio of the cross section changes with a d . As expected, more energy is transmitted for the smallest scatterers, but the transmission drops to zero as the ratio a d increases to 0.5, at which point there is no longer any separation between the cylinders.

Conclusion
In this paper, we have presented a new method of calculating the two-dimensional free-space quasi-periodic Green's function. We began by formulating the boundary integral equation for the scattering of an acoustic plane wave from a one-dimensional array of identical cylinders of arbitrary cross-section, and derived the form of the quasi-periodic Green's function. After truncating the infinite sum present in the quasi-periodic Green's function, we then used an asymptotic expansion for the remaining tail ends to produce the correction terms. This approach provides a fast and efficient way to approximate the function in a manner which is straightforward to implement. We compared our new method of calculation for the quasi-periodic Green's function against state-of-the-art methods, including Ewald's method and Kummer's transformation. We found that our results were in good agreement with the results of [11], and although less efficient than Ewald's method, the asymptotic correction is straightforward to implement and does not require the fine-tuning of multiple parameters for different frequencies. The approximation also performs well in providing results accurate up to 5 significant figures with low values of the truncation point M. In order to demonstrate the efficacy of our asymptotic correction, we implemented it as part of a BEM scheme which we have used to calculate reflection and transmission coefficients from arrays of cylinders with various cross-sections. Comparisons with previous results for circular cylinders and results generated using FEM in the COMSOL software package for the remaining geometries show good agreement, which further confirms the validity of our correction term. We have also investigated the behaviour of the transmitted energy from an array of elliptical cylinders as the size and aspect ratio of the elliptic cross-section varies. An immediate extension of the work would be to use the results of this paper to derive an asymptotic correction for the quasi-periodic Green's function that arises in the elastodynamic problem of an elastic plane wave scattering from an infinite array of cylinders. It may also be possible to extend the approach to the 3D Green's function and the scattering of a wave from a fully two-dimensional array of scatterers.