# A Direct Method for Bloch Wave Excitation by Scattering at the Edge of a Lattice. Part I: Point Scatterer Problem

A Direct Method for Bloch Wave Excitation by Scattering at the Edge of a Lattice. Part I: Point... Summary A new method for determining the reflection and transmission properties of lattices is developed. The method uses multipole expansions, and certain transformations of the algebraic equation systems that appear when boundary conditions are applied. It is more direct, and much simpler, than earlier approaches based on integral transforms and the Wiener–Hopf technique. The method is demonstrated for the case of a semi-infinite lattice of sound soft acoustic point scatterers, but can easily be generalised to account for finite size effects, and more general boundary conditions. 1. Introduction Bloch waves propagate through periodic structures without loss of energy, and have been found to exist in a number of physical contexts including electromagnetism, acoustics and elasticity (1, 2, 3). Generally, a periodic structure that supports Bloch wave propagation can only do so at certain frequencies, and in certain directions. This band structure gives rise to a range of applications, including basic components such as filters and mirrors, and more radical ideas such as optical circuits (4). Much of the literature in this area is devoted to the propagation problem, that is determining the parameter regimes in which Bloch waves can exist in a particular medium. In contrast, very little attention has been paid to the excitation problem, in which a wave incident from free space strikes the edge of a periodic structure, and some of the energy may be converted into Bloch waves. The excitation problem is generally much more difficult than the associated propagation problem, but it has enormous implications for the effectiveness of periodic media in manipulating wave propagation. For example, a band gap filter can prevent the transmission of signals at certain frequencies, but part of the signal at the desired frequencies may also be reflected back. By solving the excitation problem, it is possible to determine the proportion of incident energy that is lost in this way. A qualitative discussion of electromagnetic Bloch wave excitation is given in (4, pp. 221–227), but here the authors are only concerned with the number of diffracted Bloch waves and the directions in which these propagate. No method for determining their amplitudes is given. Perhaps the simplest case in which Bloch wave excitation can occur involves a plane wave incident on a semi-infinite lattice of acoustic point scatterers. This was considered in (5). In a subsequent paper, the same authors went on to introduce finite size effects, thereby allowing for a greater range of boundary conditions and frequencies (6). Another recent paper considers excitation of Bloch waves in a semi-infinite lattice of pinned points on a thin plate modelled by Kirchhoff theory (7). In all these three papers, the boundary value problem is converted into an infinite linear system of algebraic equations (using the multipole method (8, Chapter 4), or a reduced form for point scatterers), but these systems can only be solved directly in cases where no Bloch waves exist; otherwise truncation introduces spurious reflection effects which change the nature of the problem. However, applying an appropriate integral transform leads to a Wiener–Hopf equation. Once this is solved, the reflection and transmission properties of the lattice can be calculated. The idea of using the Wiener–Hopf method in this way was originally introduced in a study of wave diffraction by a linear array (9, 10) (see also (11)). For problems involving point scatterers, the Wiener–Hopf equation is scalar, and can be solved by a standard procedure (12). On the other hand, allowing for finite size effects leads to a matrix Wiener–Hopf equation. Although solutions to equations of this type are known to exist, obtaining them is notoriously difficult, and a general method has yet to be found (see (13) for a discussion of this issue and relevant references). In particular, the matrix equation obtained in (6) is not amenable to any of the methods recently surveyed in (14). The authors of (6) overcame this issue by approximating the field between consecutive rows as a finite expansion in grating modes. This eventually reduces the problem of solving the Wiener–Hopf equation to a complex root finding problem. However, implementing a solution to this turns out to be rather difficult, because the roots must be located with great precision, and the complexity of the function and the number of roots become greater as the frequency is increased. The purpose of this article is to introduce a more direct method for calculating the reflection and transmission properties of lattices, based on the filtering technique introduced in (15). The idea behind this is very straightforward. First, the dispersion relation from the associated propagation problem is solved, to determine the Bloch vectors for waves that may be excited inside the lattice. The algebraic equations obtained from the multipole method are then transformed in such a way that slowly convergent terms, present due to the excitation of Bloch waves, are removed. Only Bloch vectors corresponding to waves that transport energy towards infinity are included in this process; this is the radiation condition for Bloch wave excitation problems. The resulting rapidly convergent system can be truncated without introducing spurious reflection or damping effects. After inverting the truncated system numerically, the solution to the original system can be reconstructed, and all important quantities such as the Bloch wave amplitudes, and the proportions of energy transmitted into and reflected back from the lattice, can be calculated. Aside from being simpler than the Wiener–Hopf approach discussed above, our new method has the advantage that allowing for lattice elements of finite size has very little effect aside from introducing more complicated algebra. Therefore it makes sense to present the method initially for the case of acoustic point scatterers. The structure of the article is as follows. Sections 2–4 contain a review of the basic material needed for modelling wave interactions with lattices. Section 5 contains some formulas for the evaluation of quasiperiodic Green’s functions, which are required by our method. The same Green’s functions are used when the problem is solved via the Wiener–Hopf technique (see (5)), so the need for these is not a new complication. Section 6 presents a simple method for solving the propagation problem, using the quasiperiodic Green’s functions from section 5. The algebraic system of equations that governs scattering by a semi-infinite lattice is obtained in section 7. This is the same equation that was solved via the Wiener–Hopf method in (5). Indeed, all of the material in sections 2–7 can be found in the references given. Its inclusion here is necessary in order to present our method in a coherent fashion, using consistent notation. The application of the filtering method to the infinite algebraic system is performed in section 8. This is the key part of the article, and the process for transforming the system in the presence of an arbitrary number of Bloch waves is the main result. To maintain the flow of the material, a technical part of the proof that establishes the validity of the transformation is deferred to the appendix. Section 9 contains details of another idea from (15), namely infinite array subtraction. This can be used to recalculate certain coefficients in the problem, given the amplitudes of the excited Bloch waves obtained from the filtering method, and provides a stringent means of testing the accuracy of our results. Formulae for determining the reflected and transmitted fields, and a related conservation of energy criterion are obtained in sections 10 and 11. A small number of numerical results are presented in section 12. The purpose of these is to demonstrate that the method works, and to illustrate its accuracy. A greater range of results will be presented in the second part of this work, where finite size effects will be introduced. Some concluding remarks, including a detailed comparison of the filtering method and the Wiener–Hopf approach, are given in section 13. 2. Scattering by a small cylinders We will consider the interaction of two-dimensional wave fields $$u(\mathbf{r})$$ with cylinders of radius $$a$$ and infinite length (so that there is no variation in the $$z$$ direction). Outside the cylinders, the field $$u(\mathbf{r})$$ must satisfy the Helmholtz equation   $$\label{helmholtz} (\nabla^2+k^2)u(\mathbf{r}) = 0,$$ (1) and on their surfaces we have the Dirichlet boundary condition   $$\label{dirichlet} u(\mathbf{r}) = 0.$$ (2) After $$u(\mathbf{r})$$ has been determined, a time-harmonic acoustic potential can be retrieved by writing   $$U(\mathbf{r},t) = \mathrm{Re}[ u(\mathbf{r}){\mathrm e}^{-{\mathrm i}\omega t}],$$ (3) where $$\omega=ck$$, with $$\omega$$ representing frequency, and $$c$$ the speed of sound in the exterior medium. Note that $$U(\mathbf{r},t)$$ can also be interpreted as the $$z$$ component of a transverse magnetic (TM) electric field vector incident on a lattice of perfectly conducting cylinders. In this case $$c$$ represents the speed of light. Now any field incident on a cylinder centred at $$\mathbf{r}=\mathbf{R}$$ can be expanded in the form   $$u^{\mathrm i}(\mathbf{r}) = \sum_{n=-\infty}^\infty I_n {J}_n(k|\mathbf{r}-\mathbf{R}|){\mathrm e}^{{\mathrm i} n\mu},$$ (4) where $${J}_n$$ is the Bessel function of order $$n$$, and $$\mu$$ is the anticlockwise angle between the positive $$x$$ axis and the vector $$\mathbf{r}-\mathbf{R}$$. Setting $$\mathbf{r}=\mathbf{R}$$ reveals that $$I_0=u^{\mathrm i}(\mathbf{R})$$, and hence, on the surface of the cylinder we have   $$\label{expand_inc} u^{\mathrm i}(\mathbf{r}) = u^{\mathrm i}(\mathbf{R}) {J}_0(ka) + O(ka).$$ (5) If we take the scattered field to be   $$\label{single_scattered} u^{{\mathrm s}}(\mathbf{r}) = A {\mathrm{H}}_0^{(1)}(k|\mathbf{r}-\mathbf{R}|),$$ (6) where $${\mathrm{H}}_0^{(1)}$$ is a Hankel function of the first kind, then the Dirichlet boundary condition (2) is satisfied at the leading-order if the amplitude coefficient $$A$$ is given by   $$\label{foldy} A = -Z_0 u^{\mathrm i}(\mathbf{R}),$$ (7) in which   $$\label{def:Z0} Z_0 = \frac{{J}_0(ka)}{{\mathrm{H}}_0^{(1)}(ka)} = \biggl( 1 + \frac{2{\mathrm i}}{\pi}\bigl[ C + \ln(ka) - \ln 2\bigr]\biggr)^{-1} + O\bigl((ka)^2\bigr).$$ (8) Here, $$C=0.5772\ldots$$ is Euler’s constant. This method for treating scattering in the low-frequency limit dates back to (16). The same results can also be obtained using asymptotic expansions (see (8, Chapter 8)). In this way one can also derive approximate scattering coefficients for arbitrary small shapes. Higher-order terms can be obtained using shift operators (8, Sections 2.4–2.5), leading to Graf’s addition theorem and the full linear theory for multiple scattering. However, the leading-order approximation, in which $$O(ka)$$ and smaller terms are discarded, is sufficient for our purposes. 3. Two-dimensional lattices A two-dimensional lattice can be formed by placing identical objects at points with position vectors   $$\label{def:Rjp} \mathbf{R}_{jp} = j\mathbf{s}_1 + p\mathbf{s}_2, \quad j,p\in\mathbb{Z}.$$ (9) Without loss of generality, we may assume that   $$\label{def:s12} \mathbf{s}_1 = s_1 \hat{\mathbf{x}} \quad\text{and}\quad \mathbf{s}_2 = \eta_1 \hat{\mathbf{x}} + \eta_2 \hat{\mathbf{y}},$$ (10) where the circumflex denotes a unit vector,   $$s_1/2 \ge \eta_1 \ge 0 \quad\text{and}\quad \eta_2>0.$$ (11) Here we have introduced the convention that $$|\mathbf{v}|=v$$ for any vector $$\mathbf{v}$$, which will be used throughout. We will also require the reciprocal lattice vectors   $$\label{def:recip_vecs} \mathbf{s}_1^* = \frac{1}{s_1} \left(\hat{\mathbf{x}}-\frac{\eta_1}{\eta_2}\,\hat{\mathbf{y}}\right) \quad\text{and}\quad \mathbf{s}_2^* = \frac{\hat{\mathbf{y}}}{\eta_2}.$$ (12) These have the property that   $$\label{recip_prop} \mathbf{s}_j \cdot \mathbf{s}_p^* = \delta_{jp},$$ (13) where $$\delta_{jp}$$ is the Kronecker delta. Consequently, if   $$\label{def:Rmn} \mathbf{R}_{mn}^* = 2\pi(m \mathbf{s}_1^* + n \mathbf{s}_2^*), \quad m,n\in\mathbb{Z},$$ (14) then   $$\label{recip_prop2} {\mathrm e}^{{\mathrm i} \mathbf{R}_{jp}\cdot \mathbf{R}_{mn}^*} = 1,$$ (15) for all combinations of integers $$j$$, $$p$$, $$m$$ and $$n$$. Further details of lattice theory can be found in (17) and (4, Appendix B). 4. Grating modes In what follows, we will consider wave fields that have the one-dimensional quasiperiodicity property   $$\label{qp_1d} u(\mathbf{r} + s_1\hat{\mathbf{x}} ) = {\mathrm e}^{{\mathrm i} s_1 \beta_x} u(\mathbf{r}).$$ (16) A field that satisfies (16) and the Helmholtz equation (1) may be expanded as a series of grating modes. That is,   $$\label{grating_expansion} u(\mathbf{r}) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} \beta_{xj} x} \bigl[ c_j^- {\mathrm e}^{\gamma(\beta_{xj})y} + c_j^+ {\mathrm e}^{-\gamma(\beta_{xj})y} \bigr],$$ (17) where   $$\label{def:bxj} \beta_{xj} =\beta_x+2j\pi/s_1, \quad j\in\mathbb{Z},$$ (18) and   $$\label{def:gamma} \gamma(t) = \begin{cases} \sqrt{ t^2 - k^2 } &\text{if } |t|\ge k, \\ -{\mathrm i} \sqrt{ k^2 - t^2 } &\text{if } |t|< k. \end{cases}$$ (19) In general, the amplitude coefficients $$c_j^\pm$$ are different in different regions, for example, above and below a grating, or between different rows of a lattice. There are a finite number terms in (17) for which the real part of $$\gamma(\beta_{xj})$$ is zero; these represent propagating modes. Since the distinction between propagating modes and evanescent modes is important, we introduce the sets   $$\label{def:MN} \mathcal{M} = \bigl\{ j: \ j\in\mathbb{Z}, \ |\beta_{xj}|\le k \} \quad\text{and}\quad \mathcal{N} = \bigl\{ j: \ j\in\mathbb{Z}, \ |\beta_{xj}| > k \}.$$ (20) We will refer to modes with amplitude coefficient $$c_j^+$$ as upward oriented (propagating in the positive $$y$$ direction, or decaying as $$y\to\infty$$), whereas those with coefficients $$c_j^-$$ are downward oriented. To avoid certain technicalities associated with Wood’s anomalies (also sometimes called resonances), we will assume that the parameters $$\beta_x$$, $$s_1$$ and $$k$$ are such that $$|\beta_{xj}|\neq k$$ for any integer $$j$$, so that $$\gamma(\beta_{xj})\neq0$$. For the case of a linear array, a method for dealing with Wood’s anomalies was developed in (18). It is possible to adapt the analysis of semi-infinite lattices in a similar way, but this is not pertinent to the current article. For later purposes, it is useful to determine the direction in which a field of the form (17) transports energy across lines parallel to the $$x$$ axis. Now the average energy flux over any line $$\mathcal S$$ in one time-period is given by the line integral (5)   $$\label{energy_int} \langle E_{\mathcal S} \rangle = -\,\frac{P_0\omega}{2} \mathrm{Im} \int_{\mathcal S} u(\mathbf{r}) \frac{ \partial }{\partial n} u^*(\mathbf{r}) \, {\mathrm d} s,$$ (21) where the derivative is taken in a direction normal to $$\mathcal S$$ and the superscript ‘$$*$$’ denotes a complex conjugate. If $$\langle E_{\mathcal S} \rangle > 0$$, the direction of net energy flux is the same as the orientation of the normal. The constant $$P_0$$ represents the quiescent fluid density in the case of an acoustic wavefield. For a TM polarised electric field, this should be replaced by $$1/(\omega^2\mu_0)$$, where $$\mu_0$$ is the magnetic permeability of the exterior medium. Taking $$\mathcal S$$ to be the straight line from the point $$(x_0 - s_1/2,y_0)$$ to $$(x_0 + s_1/2,y_0)$$, and directing the normal parallel to $$\hat{\mathbf{y}}$$, (21) becomes   $$\label{energy2} \langle E_{\mathcal S} \rangle = -\,\frac{P_0\omega}{2} \mathrm{Im} \int_{x_0-s_1/2}^{x_0+s_1/2} u(\mathbf{r}) \frac{ \partial }{\partial y} u^*(\mathbf{r})\biggr|_{y=y_0} \, {\mathrm d} x.$$ (22) Now   $$\label{grating_int} \int_{x_0-s_1/2}^{x_0+s_1/2} {\mathrm e}^{{\mathrm i} \beta_{xj}x} {\mathrm e}^{-{\mathrm i} \beta_{xp}x} \, {\mathrm d} x = s_1 \delta_{jp},$$ (23) so if we use (17) in (22), we obtain   \begin{multline} \langle E_{\mathcal S} \rangle = -\,\frac{P_0\omega s_1}{2} \mathrm{Im} \sum_{j=-\infty}^\infty \gamma^*(\beta_{xj}) \bigl[ c_j^- {\mathrm e}^{\gamma(\beta_{xj})y_0} + c_j^+ {\mathrm e}^{-\gamma(\beta_{xj})y_0} \bigr] \\ \times \bigl[ (c_j^-)^* {\mathrm e}^{\gamma^*(\beta_{xj})y_0} - (c_j^+)^* {\mathrm e}^{-\gamma(\beta_{xj})^*y_0} \bigr]. \end{multline} (24) Considerable simplifications now occur on separating terms in which $$\gamma(\beta_{xj})$$ is imaginary from those in which it is real. We find that   $$\label{energy_prop} \langle E_{\mathcal S} \rangle = \frac{P_0\omega s_1}{2} \sum_{j\in\mathcal M} |\gamma(\beta_{xj})| \bigl( |c_j^+|^2 - |c_j^-|^2 \bigr) - P_0\omega s_1 \sum_{j\in\mathcal N}\gamma(\beta_{xj})\mathrm{Im}\bigl[ c_j^+ (c_j^-)^*\bigr],$$ (25) where the sets $$\mathcal{M}$$ and $$\mathcal{N}$$ are defined in (20). The same formula is given in (6) and also earlier (with different notation) in (19). As noted in (19), the second term on the right-hand side of (25) plays no role in problems involving a single grating illuminated by a plane wave, because in such cases there are no incoming evanescent modes, so for $$j\in\mathcal N$$ we have $$c_j^-=0$$ above the grating and $$c_j^+=0$$ below. However, this term is of crucial importance in lattice problems because there is a full set of upward- and downward-oriented grating modes between each consecutive pair of rows. Since $$x_0$$ and $$y_0$$ do not appear in (25), there is some freedom in setting the location of $$\mathcal S$$. Convenient choices are typically $$x_0=q\eta_1$$ and $$y_0=(q+1/2)\eta_2$$ for an integer $$q$$, so that there are no singularities on the integration path (see, for example, Fig. 7 in (5)). 5. Quasiperiodic Green’s functions The quasiperiodic Green’s function for a row of sources located at the points $$\mathbf{r}=\mathbf{R}_{j0}$$ (9) is given by   $$\label{def:QPG1d} G_0(\mathbf{r};\beta_x) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} js_1 \beta_x} {\mathrm{H}}_0^{(1)}(k|\mathbf{r}-\mathbf{R}_{j0}|).$$ (26) By a straightforward application of Poisson summation (20, 5), this may be rewritten in the spectral form   $$\label{spectral:QPG1d} G_0(\mathbf{r};\beta_x) = -\,\frac{2{\mathrm i}}{s_1}\sum_{j=-\infty}^\infty \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x} }{\gamma(\beta_{xj})}\, {\mathrm e}^{- \gamma(\beta_{xj}) |y| },$$ (27) where $$\beta_{xj}$$ and $$\gamma$$ are defined in (18) and (19), respectively. Closely related to $$G_0$$ is the one-dimensional lattice sum   \begin{align} \sigma_0(\beta_x) &= \lim_{r\to 0} \bigl[ G_0(\mathbf{r};\beta_x) - {\mathrm{H}}_0^{(1)}(kr) \bigr] \\ \end{align} (28)  \begin{align} &= 2\sum_{j=1}^\infty \cos(js_1\beta_x){\mathrm{H}}_0^{(1)}(js_1), \label{def:sigma0} \end{align} (29) which is sometimes called a Schlömilch series, and can be evaluated using the methods in (21). We will also need the quasiperiodic Green’s functions for multiple rows, which were originally considered in (5). For rows $$q_0,\ldots,q_1$$, we have   \begin{align} G_0^{(q_0,q_1)}(\mathbf{r};{\boldsymbol{\beta}}) &= \sum_{q=q_0}^{q_1} \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} \mathbf{R}_{jq} \cdot {\boldsymbol{\beta}} } {\mathrm{H}}_0^{(1)}(k|\mathbf{r}-\mathbf{R}_{jq}|) \label{def:22QPG} \\ \end{align} (30)  \begin{align} &= \sum_{q=q_0}^{q_1} {\mathrm e}^{{\mathrm i} q \mathbf{s}_2\cdot{\boldsymbol{\beta}}} G_0(\mathbf{r}-q\mathbf{s}_2;\beta_x), \label{22QPG2} \end{align} (31) where $${\boldsymbol{\beta}}=\beta_x\hat{\mathbf{x}}+\beta_y\hat{\mathbf{y}}$$. If $$q_0$$ and $$q_1$$ are finite integers, this has the spectral form   \begin{align} G_0^{(q_0,q_1)}(\mathbf{r};{\boldsymbol{\beta}})= -\,\frac{2{\mathrm i}}{s_1}\sum_{j=-\infty}^\infty \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x \mp \gamma(\beta_{xj})y } }{\gamma(\beta_{xj})}\,\frac{ {\mathrm e}^{q_0 w_j^{^\pm}} - {\mathrm e}^{(1+q_1)w_j^{^\pm}} } { 1 - {\mathrm e}^{w_j^{^\pm}} }, \label{MRQPG} \end{align} (32) where   \begin{align} w^{\pm}_j = \pm\eta_2 \gamma(\beta_{xj}) + {\mathrm i}( \eta_2\beta_y - 2j\pi\eta_1/s_1), \label{def:wj} \end{align} (33) and the upper and lower signs are to be taken when $$y\ge q_1\eta_2$$ and $$y\le q_0\eta_2$$, respectively. The situation in which $$q_0\eta_2 < y < q_1\eta_2$$ can be accounted for by considering the rows above and below the observer separately and using (32) twice. Finally, for semi-infinite lattices, (5) gives the formulae   \begin{align} G_0^{(q_0,\infty)}(\mathbf{r};{\boldsymbol{\beta}}) = \frac{2{\mathrm i}}{s_1} \sum_{j=-\infty}^\infty {\mathrm e}^{q_0 w_j^{^-}} \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x + \gamma(\beta_{xj})y } }{\gamma(\beta_{xj}) \big({\mathrm e}^{w_j^{^-}}-1 \big) }, \quad y \le q_0\eta_2 \label{QPG:upper} \end{align} (34) and   \begin{align} G_0^{(-\infty,q_1)}(\mathbf{r};{\boldsymbol{\beta}}) = \frac{2{\mathrm i}}{s_1} \sum_{j=-\infty}^\infty {\mathrm e}^{q_1 w_j^{^+}} \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x -\gamma(\beta_{xj})y } }{\gamma(\beta_{xj}) \big({\mathrm e}^{-w_j^{^+}}-1\big) }, \quad y \ge q_1\eta_2. \label{QPG:lower} \end{align} (35) 6. The propagation problem The problem of Bloch wave propagation through two-dimensional lattices of cylindrical Dirichlet scatterers was considered in detail in (1). Here we require only a leading-order approximation, valid when the scatterer radius is much smaller than the incident wavelength. In this case a Bloch wave has the form   $$\label{Bloch1} u({\mathbf{{r}}}) = B G_0^{(-\infty,\infty)}({\mathbf{{r}}};{\boldsymbol{\beta}}),$$ (36) for a Bloch vector $${\boldsymbol{\beta}}$$ and an arbitrary amplitude coefficient $$B$$. This has the two-dimensional quasiperiodicity property   $$\label{QP2d} u({\mathbf{{r}}}+{\mathbf{{R}}}_{jp}) = {\mathrm e}^{{\mathrm i} {\mathbf{{R}}}_{jp}\cdot {\boldsymbol{\beta}}} u({\mathbf{{r}}}),$$ (37) which incorporates the one-dimensional property (16). The Bloch vector must be chosen so that the boundary condition on the cylinder centred at the origin is satisfied. Boundary conditions elsewhere then follow by quasiperiodicity. From the perspective of this cylinder, the Bloch wave may be decomposed into incoming and scattered fields $$u_0^{\mathrm i}$$ and $$u_0^{\mathrm s}$$, where the incoming field consists of the radiation from all the other cylinders. That is,   \begin{align} u^{\mathrm i}_0({\mathbf{{r}}};{\boldsymbol{\beta}}) = B\sum_{p=-\infty}^\infty \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} {\mathbf{{R}}}_{jp}\cdot {\boldsymbol{\beta}}} {H}_0^{(1)}(k|{\mathbf{{r}}}- {\mathbf{{R}}}_{jp}|) \label{def:u0i} \end{align} (38) and   \begin{align} u^{\mathrm s}_0({\mathbf{{r}}};{\boldsymbol{\beta}}) = B {H}_0^{(1)}(kr), \end{align} (39) where the prime in (38) indicates that the term in which $$R_{jp}=0$$ is to be omitted from the summation. Comparing these to (6) and (7), we find that $$B=-Z_0 u_0^{\mathrm i}({\mathbf{{0}}};{\boldsymbol{\beta}})$$, where $$Z_0$$ is given by (8). Evaluating (38) at the origin then yields   $$\label{bloch_dispersion} 1 + Z_0 \Xi_0({\boldsymbol{\beta}}) = 0,$$ (40) where   $$\label{def:2dsigma} \Xi_0({\boldsymbol{\beta}}) = \sum_{p=-\infty}^\infty \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} {\mathbf{{R}}}_{jp}\cdot {\boldsymbol{\beta}}} {H}_0^{(1)}(k R_{jp}).$$ (41) A number of methods for computing two-dimensional lattice sums of this type are discussed in (21). One simple approach is to use the decomposition   $$\label{2dsigma_decomp} \Xi_0({\boldsymbol{\beta}}) = G_0^{(-\infty,-1)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + \sigma_0(\beta_x) + G_0^{(1,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}),$$ (42) where the terms on the right-hand side are defined in section 5. Equation (40) is the dispersion relation for Bloch waves. Given fixed values for all parameters except $$\beta_y$$, its real roots must be determined before we can proceed to consider Bloch wave excitation. To simplify this process, we use results from (21, Section 3) to conclude that $${\text{Re}}[\Xi_0]=-1$$, and so (40) becomes   $$\label{d2} 1 + W_0\,{\text{Im}}[\Xi_0({\boldsymbol{\beta}})] = 0,$$ (43) where   $$W_0 = {\mathrm i} Z_0/(1-Z_0) = {J}_0(ka) / {Y}_0(ka).$$ (44) Here, $${Y}_0$$ is the Neumann function of order zero, and the coefficient $$W_0$$ is clearly real. Now, from (9) and (10), we have $${\mathbf{{R}}}_{jp}\cdot\hat{{\mathbf{{y}}}}=p\eta_2$$, so (37) shows that all possibilities are accounted for if we consider   $$\label{beta_y_range} 0 \le \beta_y <2\pi/\eta_2.$$ (45) We must also take into account the fact that the function $$\Xi_0$$ has poles at points where $$|{\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*| = k$$, for some $$m,n\in\mathbb{Z}$$ (21, (3.18)). From (14) and (12),   $$\hat{{\mathbf{{x}}}} \cdot ({\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*) = \beta_x + \frac{2m\pi}{s_1}$$ (46) so at any pole we must have   $$\label{m_range} \left|\beta_x + \frac{2m\pi}{s_1} \right| \le k.$$ (47) Given a value for $$m$$ which satisfies (47), there are two possible values for the $$y$$ component. These are given by   $$\hat{{\mathbf{{y}}}}\cdot ({\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*) = \pm\sqrt{k^2- \bigl(\hat{{\mathbf{{x}}}}\cdot ({\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*)\bigr)^2 },$$ (48) from which we easily obtain   $$\beta_y + \frac{2n\pi}{\eta_2} = \frac{2\pi m\eta_1}{s_1 \eta_2} \pm \sqrt{k^2 - (\beta_x+2m\pi/s_1)^2 }.$$ (49) Since $$n\in\mathbb{Z}$$ it now follows that, for each $$m$$ satisfying (47), there are precisely two poles such that $$\beta_y$$ satisfies (45), unless the square root evaluates to zero, in which case there is only one pole for this particular choice of $$m$$. Another useful piece of information that can be obtained from the propagation problem is the direction in which each Bloch wave carries energy across lines parallel to the $$x$$ axis. Now Bloch waves do not have well-defined phase velocities, because integer multiples of $${\mathbf{{s}}}_1^*$$ and $${\mathbf{{s}}}_2^*$$ can be added to $${\boldsymbol{\beta}}$$ with no substantive effect in view of (15) and (37). Therefore, we must use information about energy transportation to determine the direction of propagation (see (4, pp. 222–225) for further discussion). To this end, we evaluate (25) for the Bloch wave represented by (36). Thus, using (34) and (35) with $$q_1=q$$, and $$q_0=1+q$$ we find that for $$q \eta_2 \le y \le (1+q) \eta_2$$, a Bloch wave can be expanded in the form (17), with   $$c_j^+ = \frac{2{\mathrm i} B {\mathrm e}^{q w_j^+}}{s_1\gamma(\beta_{xj})( {\mathrm e}^{-w_j^+} - 1 )} \quad\text{and}\quad c_j^- = \frac{2{\mathrm i} B {\mathrm e}^{(1+q) w_j^-}}{s_1\gamma(\beta_{xj})( {\mathrm e}^{w_j^-} - 1 )}.$$ (50) If $$\gamma(\beta_{xj})$$ is imaginary, then so too are $$w_j^+$$ and $$w_j^-$$, and so   $$\label{cjsq} |c_j^+|^2 = \frac{4 |B|^2}{s_1^2\bigl|\gamma(\beta_{xj})( {\mathrm e}^{-w_j^+} - 1 )\bigr|^2} \quad\text{and}\quad |c_j^-|^2 = \frac{4 |B|^2}{s_1^2\bigl|\gamma(\beta_{xj})( {\mathrm e}^{w_j^-} - 1 )\bigr|^2}, \quad j \in\mathcal{M}.$$ (51) On the other hand, if $$\gamma(\beta_{xj})$$ is real, then $$(w_j^-)^* = -w_j^+$$, so   $$\label{cjpcjm} {\text{Im}}\bigl[ c_j^+ (c_j^-)^*\bigr] = \left|\frac{B}{s_1 \gamma(\beta_{xj})}\right|^2 {\text{Im}} \bigl[{\text{csch}}^2( w_j^+/2)\bigr], \quad j \in\mathcal{N}.$$ (52) There is no dependence on $$q$$ here, which should be expected because Bloch waves propagate through the lattice without loss of energy. For a given value for $$\beta_y$$, we may now compute $$\langle E_{\mathcal S}\rangle$$ using (25), up to a factor of $$|B|^2$$, which plays no role in determining the direction of propagation. If $$\langle E_{\mathcal S} \rangle > 0$$, or $$\langle E_{\mathcal S} \rangle < 0$$, then the Bloch wave transports energy in the direction of increasing, or decreasing $$y$$, respectively. A Bloch wave with $$\langle E_{\mathcal S} \rangle=0$$ does not transport any energy in the $$y$$ direction. 7. The excitation problem We now consider a plane wave   $$u^{\mathrm i}({\mathbf{{r}}}) = {\mathrm e}^{{\mathrm i} k(x\cos\psi_0 +y\sin\psi_0)}$$ (53) incident on a semi-infinite lattice formed from small cylinders centred at the points $${\mathbf{{r}}}={\mathbf{{R}}}_{jp}$$, with $$j\in\mathbb{Z}$$ and $$p=0,1,\ldots$$ (see Fig. 1). The scattered field for this problem inherits the one-dimensional quasiperiodicity property (16) with $$\beta_x=k\cos\psi_0$$, and so we have   $$\label{def:us} u^{\mathrm s}({\mathbf{{r}}}) = \sum_{p=0}^\infty A_p \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} j s_1k\cos\psi_0} {H}_0^{(1)}(k|{\mathbf{{r}}}-{\mathbf{{R}}}_{jp}|),$$ (54) for as yet unknown coefficients $$A_p$$. It should be noted that the quasiperiodicity condition generally precludes the existence of surface waves along the edge at $$y=0$$, because $$0\in\mathcal{M}$$ for all angles of incidence (see (20)), meaning there is always at least one grating mode that transports energy into free space. Fig. 1. View largeDownload slide A plane wave incident on a semi-infinite lattice Fig. 1. View largeDownload slide A plane wave incident on a semi-infinite lattice We will apply the Dirichlet boundary condition (2) on the cylinders centred at the points $${\mathbf{{r}}}={\mathbf{{R}}}_{0q}$$, for $$q=0,1,\ldots$$ Boundary conditions elsewhere then follow by quasiperiodicity. The field incoming towards the cylinder located at $${\mathbf{{r}}}={\mathbf{{R}}}_{0q}$$, which consists of the incident wave and the radiation from all the other cylinders, is given by   $$\label{uiq} u^{\mathrm i}_q({\mathbf{{r}}}) = {\mathrm e}^{{\mathrm i} k(x\cos\psi_0 +y\sin\psi_0)} + \sum_{p=0}^\infty A_p \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} j s_1 k\cos\psi_0} {H}_0^{(1)}(k|{\mathbf{{r}}}-{\mathbf{{R}}}_{jp}|),$$ (55) where the prime indicates that the term in which $$j=0$$ and $$p=q$$ is omitted from the summation. The scattered response is   $$u^{\mathrm s}_q({\mathbf{{r}}}) = A_q{H}_0^{(1)}(k|{\mathbf{{r}}}-{\mathbf{{R}}}_{0q}|).$$ (56) Comparing these to (6) and (7), we find that $$A_q=-Z_0 u_q^{\mathrm i}({\mathbf{{r}}}_{0q})$$. Evaluating (55) at $${\mathbf{{r}}}={\mathbf{{r}}}_{0q}$$ then yields   $$\label{lin_sys1} A_q + Z_0 \sum_{p=0}^\infty A_p \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} j s_1 k\cos\psi_0} {H}_0^{(1)}(kR_{j,p-q}) = T_q, \quad q = 0,1,\ldots$$ (57) where   $$T_q = -Z_0 {\mathrm e}^{{\mathrm i} q k(\eta_1 \cos\psi_0 + \eta_2\sin\psi_0)}.$$ (58) If we now write   $$\label{def:Sp} S_p = \left\{\begin{array}{ll} \sigma_0(k\cos\psi_0) &\text{if } p = 0,\\[3pt] G_0(p{\mathbf{{s}}}_2;k\cos\psi_0) &\text{otherwise,} \end{array}\right.$$ (59) then (57) becomes   $$\label{lin_sys2} A_q + Z_0 \sum_{p=0}^\infty A_p S_{q-p} = T_q, \quad q = 0,1,\ldots$$ (60) In cases where no Bloch waves are excited, $$A_q\to0$$ as $$q\to\infty$$, and this system can be solved by truncation. No further conditions are needed. However, if Bloch waves are excited then $$\hat{A}_q\not\to 0$$, and (60) cannot be solved by truncation. Moreover, the solution is not unique, because we have yet to apply a radiation condition. The appropriate condition is to forbid Bloch modes that transport energy towards $$y=0$$ from appearing in the far field inside the lattice. This will be incorporated into the analysis as part of the filtering method, discussed in the next section. 8. The filtering method In a case where a single Bloch wave is excited, we can write   $$A_q = B {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}} + \hat{A}_q.$$ (61) The correct value for $$\beta_y$$ must yield a solution to (43) with $$\langle E_{\mathcal S}\rangle \ge 0$$ as described in section 6. This has two effects: it causes $$\hat{A}_q$$ to vanish as $$q\to\infty$$, and it ensures that the far field inside the lattice does not include any Bloch waves incoming from infinity, so that the radiation condition is satisfied. Now the decomposition (61) introduces a new unknown into the problem (the Bloch wave amplitude $$B$$), but this can be ‘filtered out’ by writing   $$A_q^{(1)} = \begin{cases} A_0 &\text{if } q = 0, \\ A_q - {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}}A_{q-1} &\text{otherwise}. \end{cases}$$ (62) For $$q>0$$, we then have $$A_q^{(1)} = \hat{A}_q - {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}}\hat{A}_{q-1}$$, which vanishes as $$q\to\infty$$. Solving the recurrence relation for $$A_q$$, we find that   $$A_q = \sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}}.$$ (63) There are now two ways in which we may obtain a system of equations for $$A_q^{(1)}$$. One possibility is to take phase-shifted differences between consecutive equations in (60). This was the preferred approach in (15), but the equation in which $$q=0$$ requires special treatment (because there is no preceding equation). If two Bloch waves are excited, the equations with $$q=0$$ and $$q=1$$ must be treated as special cases, and so on, meaning that the process is difficult to generalise to account for arbitrary numbers of Bloch waves. Therefore we use the alternative method (also suggested in (15)), which is to substitute (63) into (60). In this way, we obtain   $$\sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} + Z_0 \sum_{p=0}^\infty \sum_{j=0}^p A_j^{(1)} {\mathrm e}^{{\mathrm i} (p-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} S_{q-p} = T_q, \quad q=0,1,\ldots$$ (64) Again, this cannot be solved by truncation, because the inner sum in the second term involves $$A_0^{(1)},A_1^{(1)},\ldots$$ for every $$p$$. However, if we interchange the summations, we find that   $$\sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} + Z_0 \sum_{j=0}^\infty A_j^{(1)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} \Gamma_{j-q}({\boldsymbol{\beta}}) = T_q,$$ (65) where   $$\Gamma_j({\boldsymbol{\beta}}) =\left\{ \begin{array}{ll} G_0^{(j,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}) &\text{if } j > 0,\\[3pt] G_0^{(1,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + \sigma_0(k\cos\psi_0) &\text{if } j = 0, \\[3pt] G_0^{(j,-1)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + G_0^{(1,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + \sigma_0(k\cos\psi_0) &\text{if } j < 0. \end{array}\right.$$ (66) Since $$A_q^{(1)}\to0$$ as $$q\to\infty$$, equation (65) can be truncated at $$j=P$$ and $$q=P$$, say. The resulting finite system is then solved numerically, and (63) is used to compute the original coefficients. If $$P$$ is such that the discarded coefficients are numerically negligible, then the Bloch wave amplitude can be calculated using (61); thus   $$B = {\mathrm e}^{-{\mathrm i} P {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}} A_P.$$ (67) We now look to generalise the above procedure, to account for cases in which multiple Bloch waves are excited. In place of (61), we write   $$A_q = \hat{A}_q + \sum_{m=1}^\lambda B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(m)}}.$$ (68) Again, the correct choices for the $$y$$ components of the Bloch vectors are those values in the range (45) that yield solutions to (43) with $$\langle E_{\mathcal S}\rangle \ge 0$$, so that $$\hat{A}_q\to0$$ as $$q\to\infty$$, and the far field cannot include incoming Bloch waves. We also define   $$A_q^{(\lambda)} =\left\{ \begin{array}{ll} A_q & \text{if } q=0\text{ or } \lambda = 0, \\[3pt] A_q^{(\lambda-1)} - {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(\lambda)}} A_{q-1}^{(\lambda-1)} & \text{otherwise,} \end{array}\right.$$ (69) so that $$A_q^{(0)}$$ is the original unfiltered coefficient, whereas $$A_q^{(\lambda)}$$ has the first $$\lambda$$ Bloch waves filtered out for $$q\ge\lambda$$. Solving the recurrence relation for $$A_q^{(\lambda-1)}$$ yields   $$A_q^{(\lambda-1)} = \sum_{j=0}^q A_j^{(\lambda)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda)}},$$ (70) but we require an expression for $$A_q^{(0)}$$ in terms of $$A_q^{(\lambda)}$$, for substitution into the linear system (60). If we start with $$A_q^{(0)}$$ and use (70) twice, we obtain   \begin{align*} A_q^{(0)} &= \sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} \\ &= \sum_{j=0}^q \sum_{p=0}^j A_p^{(2)} {\mathrm e}^{{\mathrm i}(j-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} \\ &= \sum_{p=0}^q A_p^{(2)} {\mathrm e}^{-{\mathrm i} p{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} \sum_{j=p}^q {\mathrm e}^{{\mathrm i} j {\mathbf{{s}}}_2\cdot({\boldsymbol{\beta}}^{(2)}-{\boldsymbol{\beta}}^{(1)})} \\ & = d_{2,1}\sum_{p=0}^q A_p^{(2)} {\mathrm e}^{{\mathrm i} (q-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} + d_{1,2}\sum_{p=0}^q A_p^{(2)} {\mathrm e}^{{\mathrm i}(q-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}}, \end{align*} where   $$d_{n,m} = \frac{1}{1-{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot({\boldsymbol{\beta}}^{(n)}-{\boldsymbol{\beta}}^{(m)})}}.$$ (71) Repeating this process, we find that   \begin{align*} A_n^{(0)} &= d_{2,1}\sum_{p=0}^q \sum_{j=0}^p A_j^{(3)} {\mathrm e}^{{\mathrm i}(p-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(3)}} {\mathrm e}^{{\mathrm i} (q-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} + d_{1,2}\sum_{p=0}^q \sum_{j=0}^p A_j^{(3)} {\mathrm e}^{{\mathrm i}(p-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(3)}} {\mathrm e}^{{\mathrm i}(q- p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} \\ &= d_{2,1}d_{3,1} \sum_{j=0}^q A_j^{(3)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} + d_{1,2} d_{3,2}\sum_{j=0}^q A_j^{(3)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} \\ &\quad + d_{1,3}d_{2,3} \sum_{j=0}^q A_j^{(3)} {\mathrm e}^{{\mathrm i}(q- j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(3)}}, \end{align*} since   \begin{equation*} d_{2,1}d_{1,3}+d_{1,2}d_{2,3} = d_{1,3}d_{2,3}. \end{equation*} Therefore we may postulate the general result   $$A_q^{(0)} = \sum_{m=1}^\lambda Q_m^\lambda \sum_{j=0}^q A_j^{(\lambda)} {\mathrm e}^{{\mathrm i} (q-j) {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}},$$ (72) where   $$Q_1^1 = 1 \quad\text{and}\quad Q_m^\lambda = \prod_{\substack{n=1\\ n\neq m}}^\lambda d_{n,m}, \quad \lambda > 1.$$ (73) This is correct in the cases where $$\lambda=1$$, $$2$$ and $$3$$. Let us attempt an inductive step by using (70) in (72). We find that   \begin{align*} A_q^{(0)} &= \sum_{m=1}^\lambda Q_m^\lambda \sum_{j=0}^q \sum_{p=0}^j A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(j-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}} {\mathrm e}^{{\mathrm i} (q-j) {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}} \\ &= \sum_{m=1}^\lambda Q_m^\lambda {\mathrm e}^{{\mathrm i} q {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} \sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{-{\mathrm i} p {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}} \sum_{j=p}^q {\mathrm e}^{{\mathrm i} j{\mathbf{{s}}}_2\cdot({\boldsymbol{\beta}}^{(\lambda+1)}-{\boldsymbol{\beta}}^{(m)})} \\ &= \sum_{m=1}^\lambda Q_m^\lambda d_{\lambda+1,m}\sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} + \sum_{m=1}^\lambda Q_m^\lambda d_{m,\lambda+1} \sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}} \\ &= \sum_{m=1}^\lambda Q_m^{\lambda+1}\sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} + \biggl( \sum_{m=1}^\lambda Q_m^\lambda d_{m,\lambda+1} \biggr)\sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}}. \end{align*} If we can show that   $$\sum_{m=1}^\lambda Q_m^\lambda d_{m,\lambda+1} = Q_{\lambda+1}^{\lambda+1},$$ (74) then (72) is established for all $$\lambda\in\mathbb{N}$$. Evidently (74) is true if $$\lambda=1$$. Otherwise it is equivalent to writing   $$\sum_{m=1}^\lambda \prod_{\substack{n=1 \\ n\neq m}}^\lambda \frac{ d_{n,m} }{ d_{n,\lambda+1}} = 1,$$ (75) and if we set $$z_m={\mathrm e}^{-{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}$$, the proof follows immediately on applying the result in the appendix. We may now substitute (72) into (60), recalling that $$A_q^{(0)}=A_q$$. In this way we find that the linear system for $$A_q^{(\lambda)}$$ is   $$\sum_{m=1}^\lambda Q_m^\lambda\left[ \sum_{j=0}^q A_j^{(\lambda)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(m)}} + Z_0 \sum_{j=0}^\infty A_j^{(\lambda)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(m)}} \Gamma_{j-q}({\boldsymbol{\beta}}^{(m)})\right] = T_q, \quad q=0,1,\ldots$$ (76) This can be solved by truncation at $$j=P$$ and $$q=P$$, say, and the original coefficients can be computed using (72). We may also determine the Bloch wave amplitudes using (72). Thus, if $$A_j^{(\lambda)}$$ is negligible for $$j\ge P$$, then   $$A_q^{(0)} = \sum_{m=1}^\lambda \biggl( Q_m^\lambda \sum_{j=0}^P A_j^{(\lambda)} {\mathrm e}^{-{\mathrm i} j {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}}\biggr) {\mathrm e}^{{\mathrm i} q {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}}, \quad q \ge P,$$ (77) from which we can read off the Bloch wave amplitudes as   $$B^{(m)} = Q_m^\lambda \sum_{j=0}^P A_j^{(\lambda)} {\mathrm e}^{-{\mathrm i} j {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}}.$$ (78) In the case where $$\lambda=1$$ (so that $$m=1$$ is the only possibility), using (73) and (70) reduces this to   $$B^{(1)} = \sum_{j=0}^P A_j^{(\lambda)} {\mathrm e}^{-{\mathrm i} j {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(1)}}} = {\mathrm e}^{-{\mathrm i} P {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(1)}} A_P^{(0)},$$ (79) in agreement with (67). 9. Infinite array subtraction Another idea from (15) offers a useful means of testing our results. If the Bloch vectors $${\boldsymbol{\beta}}^{(m)}$$ and associated amplitude coefficients $$B^{(m)}$$ have been determined, then we may form a linear system of equations for $$\hat{A}_q$$ by substituting (68) into (60), and taking known terms to the right-hand side. In this way, we obtain   $$\hat{A}_q + Z_0 \sum_{p=0}^\infty \hat{A}_p S_{q-p} = T_q - \sum_{m=1}^\lambda B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)} } \Bigl[ 1 + Z_0 \Gamma_{-q}\bigl({\boldsymbol{\beta}}^{(m)}\bigr)\Bigr], \quad q = 0,1,\ldots$$ (80) where $$\Gamma$$ is defined in (66). In this way, the slowly decaying terms are subtracted from the linear system, leaving a rapidly convergent system for the coefficients $$\hat{A}_q$$. The new system can be simplified slightly by exploiting the fact that each Bloch wave is a solution to the propagation problem. Thus, the dispersion relation (40) can be written in the form   $$1 + Z_0 \bigl[ G_0^{(-\infty,-q-1)}(0;{\boldsymbol{\beta}}) + \Gamma_{-q}({\boldsymbol{\beta}}) \bigr] = 0,$$ (81) so that   $$\hat{A}_q + Z_0 \sum_{p=0}^\infty \hat{A}_p S_{q-p} = T_q + Z_0 \sum_{m=1}^\lambda B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} G_0^{(-\infty,-q-1)}\bigl({\mathbf{{0}}};{\boldsymbol{\beta}}^{(m)}\bigr), \quad q = 0,1,\ldots$$ (82) It should be emphasised that this system will only converge if the correct values for $${\boldsymbol{\beta}}^{(m)}$$ and $$B^{(m)}$$ are inserted on the right-hand side. 10. Reflection and transmission It is convenient to represent the reflected and transmitted fields in terms of scattering angles $$\psi_j$$, which are defined so that   $$k\cos\psi_j = k\cos\psi_0 + 2j\pi/s_1 \quad\text{and}\quad k\sin\psi_j = {\mathrm i}\gamma(k\cos\psi_j),$$ (83) where $$\gamma$$ is given by (19). We also introduce the quantities   $$\rho_j = {\mathrm e}^{-{\mathrm i} k( \eta_1\cos\psi_j+\eta_2\sin\psi_j )} \quad\text{and}\quad \tau_j = {\mathrm e}^{{\mathrm i} k(\eta_1\cos\psi_j-\eta_2 \sin\psi_j)}.$$ (84) Note that $$\cos\psi_j$$ is always real, whereas $$\sin\psi_j$$ may be positive real or positive imaginary (but not zero since we are not dealing with Wood’s anomalies; see section 4). Consequently, $$|\rho_j|=|\tau_j|\to\infty$$ as $$|j|\to\infty$$. Using this notation, we can simplify the exponentials involving $$w_j^\pm$$, which appear in the spectral representations for the quasiperiodic Green’s functions (34) and (35). Thus, since $$\beta_x=k\cos\psi_0$$ for all Bloch vectors $${\boldsymbol{\beta}}$$, we find that   $${\mathrm e}^{w_j^+} = \rho_j {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}} \quad\text{and}\quad {\mathrm e}^{w_j^-} = \tau_j^{-1} {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}}.$$ (85) Below the lattice, the total field can be expanded as a series of grating modes as in (17), with $$\beta_{xj}=k\cos\psi_j$$. The only upward-oriented mode in this region is the incident field, so   $$u({\mathbf{{r}}}) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} kx\cos\psi_j} \bigl[ c_{j0}^- {\mathrm e}^{-{\mathrm i} ky\sin\psi_j} + \delta_{j0} {\mathrm e}^{{\mathrm i} ky\sin\psi_j} \bigr], \quad y \le 0,$$ (86) where $$\delta_{j0}=1$$ if $$j=0$$ and $$0$$ otherwise. To determine $$c_{j0}^-$$, we insert the quasiperiodic Green’s function (26) into the representation of the scattered field (54), to obtain   $$u^{\mathrm s}({\mathbf{{r}}}) = \sum_{p=0}^\infty A_p G_0({\mathbf{{r}}}-p{\mathbf{{s}}}_2;k\cos\psi_0).$$ (87) Then, introducing the decomposition (68) and using (31), we find that   $$u^{\mathrm s}({\mathbf{{r}}}) = \sum_{m=1}^\lambda B^{(m)} G_0^{(0,\infty)}\bigl({\mathbf{{r}}};{\boldsymbol{\beta}}^{(m)}\bigr) + \sum_{p=0}^\infty \hat{A}_p G_0({\mathbf{{r}}}-p{\mathbf{{s}}}_2;\beta_x).$$ (88) Finally, after using the spectral representations (27) and (34) along with (83) and (85), we may read off the reflection coefficients as   $$c_{j0}^- = {}\frac{2}{ks_1 \sin\psi_j}\Biggl[ \sum_{m=1}^\lambda \frac{B^{(m)}\tau_j }{ \tau_j-{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}} + \sum_{p=0}^\infty \hat{A}_p \tau_j^{-p} \Biggr].$$ (89) Evidently, the far field below the lattice may be obtained by simply limiting the summation in (86) to $$j\in\mathcal M$$, so that $$\sin\psi_j$$ is real, meaning that propagating modes are retained, whereas evanescent modes are discarded. Calculation of the transmitted field is more complicated. Between each consecutive pair of rows, the total field can again be represented in the form (17), this time with a full set of upward- and downward-oriented modes. That is,   $$u({\mathbf{{r}}}) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} kx\cos\psi_j} \bigl[ c_{jq}^- {\mathrm e}^{-{\mathrm i} ky\sin\psi_j} + c_{jq}^+ {\mathrm e}^{{\mathrm i} ky\sin\psi_j} \bigr], \quad (q-1)\eta_2 \le y \le q\eta_2, \quad q \in\mathbb{N}.$$ (90) Since the far field inside the lattice can consist only of Bloch waves, we decompose the total field into the form   $$u({\mathbf{{r}}}) = \hat{u}({\mathbf{{r}}}) + u^{\mathrm{b}}({\mathbf{{r}}}),$$ (91) with   \begin{align} \hat{u}({\mathbf{{r}}}) = {\mathrm e}^{{\mathrm i} k(x\cos\psi_0+y\sin\psi_0)} - \sum_{m=1}^{\lambda} B^{(m)} G_0^{(-\infty,-1)}\bigl({\mathbf{{r}}};{\boldsymbol{\beta}}^{(m)}\bigr) + \sum_{p=0}^\infty \hat{A}_p G_0({\mathbf{{r}}}-p{\mathbf{{s}}}_2;\beta_x) \end{align} (92) and   \begin{align} u^{\mathrm{b}}({\mathbf{{r}}}) = \sum_{m=1}^{\lambda} B^{(m)}G_0^{(-\infty,\infty)}\bigl({\mathbf{{r}}};{\boldsymbol{\beta}}^{(m)}\bigr). \end{align} (93) We then write   $$c_{jq}^\pm = \hat{c}_{jq}^\pm + b_{jq}^\pm,$$ (94) where the terms on the right-hand side are the contributions from (92) and (93), respectively. For $$\hat{c}_{jq}^-$$, we observe that only the last term on the right-hand side of (92) includes downward-oriented modes (and only for $$p\ge q$$). Therefore, on using (27), we find that   $$\hat{c}_{jq}^- = \frac{2\tau_j^{-q}}{ks_1\sin\psi_j}\sum_{p=0}^\infty \hat{A}_{p+q} \tau_j^{-p}.$$ (95) One the other hand, all three terms in (92) contribute to $$\hat{c}_j^+$$. Again using (27), and also (35), we obtain   $$\hat{c}_{jq}^+ = \delta_{j0} + \frac{2}{ks_1\sin\psi_j}\biggl[ \sum_{m=1}^\lambda \frac{B^{(m)}} {1-\rho_j {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}} + \sum_{p=0}^{q-1} \hat{A}_p \rho_j^p \biggr].$$ (96) Similarly, for the contributions to $$u^\mathrm{b}$$, we use (34) with $$q_0=q$$ and (35) with $$q_1=q-1$$ to show that   $$b_{jq}^{+} = \frac{2\rho_j^q}{ks_1\sin\psi_j} \sum_{m=1}^\lambda \frac{B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}{\rho_j{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}-1} \quad\text{and}\quad b_{jq}^{-} = \frac{2\tau_j^{-q}}{ks_1\sin\psi_j}\sum_{m=1}^\lambda \frac{ B^{(m)}{\mathrm e}^{{\mathrm i} q {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} } {1-\tau_j^{-1}{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}.$$ (97) We now consider the various contributions to the far field inside the lattice. It should be noted that, because each instance of (90) (that is, with a particular value for $$q$$) is valid for a finite range of $$y$$ values, we cannot simply discard evanescent modes. The behaviour of the amplitude coefficients as $$q\to\infty$$ must also be taken into account. Let us begin with (95). We write   $$y = (q-1) \eta_2 + y_0, \quad 0 \le y_0 \le \eta_2,$$ (98) so that the modes have the form   $$\hat{c}_{jq}^- {\mathrm e}^{{\mathrm i} k (x\cos\psi_j-y\sin\psi_j)} = \frac{2{\mathrm e}^{-{\mathrm i} q k \eta_1\cos\psi_j} }{ks_1\sin\psi_j}\sum_{p=0}^\infty \hat{A}_{p+q} \tau_j^{-p} {\mathrm e}^{{\mathrm i} k (x\cos\psi_j+(\eta_2-y_0)\sin\psi_j)}.$$ (99) This disappears as $$q\to\infty$$, due to the decay in the coefficients $$\hat{A}_{p+q}$$. Next consider (96), in the case where $$j\in\mathcal N$$. Since $$0\in\mathcal{M}$$, the first term on the right-hand side disappears. The second term has no dependence on $$q$$, and makes no contribution to the far field because it is to be multiplied by $${\mathrm e}^{-{\mathrm i} k y \sin\psi_j}$$, which decays exponentially as $$y$$ is increased. For the third term, we can proceed as above using (98), and write   $$\sum_{p=0}^{q-1} \hat{A}_p \rho_j^p \, {\mathrm e}^{{\mathrm i} k(x\cos\psi_j+y\sin\psi_j)} = {\mathrm e}^{{\mathrm i} k(x-(q-1)\eta_1)\cos\psi_j} {\mathrm e}^{{\mathrm i} k y_0\sin\psi_j} \sum_{p=0}^{q-1} \hat{A}_p \rho_j^{p-q+1}.$$ (100) In the far field, these contributions may be small because they originate from some distance away and decay exponentially in $$y$$, or, if they originate from nearby rows, because the amplitude coefficient $$\hat{A}_p$$ multiplying the quasiperiodic Green’s function is negligible. To see this mathematically, we split the series after roughly half the required terms. Thus,   $$\sum_{p=0}^{q-1} \hat{A}_p \rho_j^{p-q+1} = \sum_{p=0}^{{[(q-1)/2]}} \hat{A}_p \rho_j^{p-q+1} + \sum_{{p=[(q-1)/2]+1}}^{q-1} \hat{A}_p \rho_j^{p-q+1},$$ (101) where $$[\cdot]$$ means integer part (round towards $$-\infty$$). The first half of the sum disappears as $$q\to\infty$$ because $$|\rho_j|>1$$, and the second half disappears because $$|\rho_j|^{p-q+1}\le 1$$ and $$\hat{A}_p\to0$$ as $$p\to\infty$$. Finally, consider (96) for propagating modes. Since the far field consists of Bloch waves only, it must be the case that   $$\delta_{j0} + \frac{2}{ks_1\sin\psi_j}\biggl[ \sum_{m=1}^\lambda \frac{B^{(m)}} {1-\rho_j {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}} + \sum_{p=0}^\infty \hat{A}_p \rho_j^p \biggr] = 0, \quad j\in\mathcal{M}.$$ (102) In (5), the fact that these terms vanish identically arises directly from the analysis. Despite considerable effort on the part of the authors, we have not been able to prove this via the current method. Instead, we must view (102) as a necessary but not sufficient condition for the correctness of our numerical results, much like the conservation of energy condition derived in the next section. Ultimately, the far field as $$y\to\infty$$ is given by (90), with $$c_{jq}^\pm$$ replaced by $$b_{jq}^\pm$$ (that is, discarding $$\hat{c}_{jq}^\pm$$ from (94)). 11. Conservation of energy A conservation of energy condition for this problem was derived in (5), and we need not repeat the details here. The procedure is to evaluate the integral (21) around the parallelogram with vertices at   $${\mathbf{r}} = {}\pm \tfrac{1}{2}{\mathbf{s}}_1 \pm (N+\tfrac{1}{2}) {\mathbf{s}}_2, \quad N\in\mathbb{N},$$ (103) taking the derivative in the direction of the outgoing normal. The contributions from the sides parallel to $${\mathbf{s}}_2$$ cancel each other, due to the one-dimensional quasiperiodicity property (16). For the edge beneath the lattice, we simply evaluate (25) for the reflected field (86), taking into account the fact that the outgoing normal is $$-\hat{{\mathbf{y}}}$$, so that a factor of $$-1$$ is introduced. Since $$0\in\mathcal M$$ for any angle of incidence, the sum over $$\mathcal{N}$$ disappears, and we find that the time averaged energy flux across this line is   $$\label{E1} E_1 = {}-\frac{P_0\omega s_1 k}{2} \biggl[ \sin\psi_0 - \sum_{j\in\mathcal M} \sin\psi_j |c_{j0}^-|^2 \biggr],$$ (104) where $$c_{j0}^-$$ is given by (89). If no Bloch waves are excited, then all of the energy incident on the lattice must be reflected back, so we must have $$E_1=0$$. If Bloch waves are excited, then the average energy flux is in direction of the inward normal (some incident wave energy is transmitted into the lattice, so less is reflected back). In this case $$E_1<0$$ and   $$\label{cons_energy} E_1 + E_2 = 0,$$ (105) where $$E_2$$ is the time averaged energy flux across the upper edge of the parallelogram. To determine $$E_2$$, we again use (25), this time with the outgoing normal directed upward. We may assume that $$N$$ is sufficiently large to allow contributions from $$\hat{c}_{jq}^+$$ to be neglected (see section 10, above), so that the amplitude coefficients are given by (97). Some simplification then occurs because $$|\rho_j|=|\tau_j|=1$$ if $$j\in\mathcal{M}$$, whereas if $$j\in\mathcal{N}$$ then $$\tau_j^*=\rho_j$$. The average energy flux across the upper edge is therefore given by   \begin{align}\label{E2} E_2 &= \frac{2P_0\omega }{k s_1} \sum_{j\in\mathcal M} \frac{1}{\sin\psi_j} \left\{ \left| \sum_{m=1}^\lambda \frac{B^{(m)} {\mathrm e}^{{\mathrm i} N{\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}{\rho_j{\mathrm e}^{{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}-1}\right|^2 - \left|\sum_{m=1}^\lambda \frac{ B^{(m)}{\mathrm e}^{{\mathrm i} N{\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}} } {1-\tau_j^{-1}{\mathrm e}^{{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}\right|^2 \right\} \nonumber\\ &\quad - \frac{4P_0\omega}{k s_1} \sum_{j\in\mathcal N}\frac{1}{|\sin\psi_j|}\Im\left\{ \sum_{m=1}^\lambda \frac{B^{(m)} {\mathrm e}^{{\mathrm i} N{\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}{\rho_j{\mathrm e}^{{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}-1} \sum_{m=1}^\lambda \frac{ (B^{(m)})^*{\mathrm e}^{-{\mathrm i} N {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}} } {1-\rho_j^{-1}{\mathrm e}^{-{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}\right\}. \end{align} (106) If a single Bloch wave is excited, then $$\lambda=1$$ and the factors involving $$N$$ all cancel. The meaning of this is important: we assumed $$N$$ to be large so that terms that do not contribute to the far field could be neglected. The resulting expression then consists only of contributions from Bloch waves. It does not depend on $$N$$ because Bloch waves propagate through the lattice without loss of energy. If multiple Bloch waves are excited then $$\lambda>1$$, and the explicit cancellation no longer occurs. However, $$E_2$$ still consists entirely of contributions from Bloch waves, and cannot depend on $$N$$. Therefore (106) can always be simplified by setting $$N=0$$, a fact which we have confirmed numerically. 12. Numerical results Figures 2 and 3 show the magnitudes of the coefficients $$A_j$$, $$A_j^{(\lambda)}$$ and $$\hat{A}_j$$ for parameters that lead to different behaviours as $$j$$ is increased. All of the results were computed using double precision arithmetic ($${\sim}16$$ decimal digits). For cases where Bloch waves are excited, a stringent test on accuracy is performed by computing $$A_j$$ using (76) and (72), and again by infinite array subtraction using (82) and (68), retaining only the Bloch wave amplitudes $$B^{(m)}$$ and vectors $${\boldsymbol{\beta}}^{(m)}$$ from the earlier calculation. We also confirm that (105) and (102) are satisfied in each case. Some care must be taken in assessing the accuracy of these, since the exact value is zero in both cases, so that the relative error is undefined. For (105), we use the number of cancelling digits in $$E_1$$ and $$E_2$$, or in cases where no Bloch waves are excited, the number of cancelling digits in the two terms in (104). For (102), we use the fact that there exists a propagating mode of unit amplitude (the incident field), so the number of correct significant figures is one greater than the number of correct decimal places. Fig. 2 View largeDownload slide Magnitudes of coefficients for a rectangular lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0,1]$$, $$a=0.005$$, $$\psi_0=0.25\pi$$ and (a) $$k=1.5$$, (b) $$k=3.0$$ Fig. 2 View largeDownload slide Magnitudes of coefficients for a rectangular lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0,1]$$, $$a=0.005$$, $$\psi_0=0.25\pi$$ and (a) $$k=1.5$$, (b) $$k=3.0$$ Fig. 3 View largeDownload slide Magnitudes of coefficients for a skewed lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0.1,1.2]$$, $$a=0.005$$, and (a) $$\psi_0=0.15\pi$$, $$k=3.7$$, (b) $$\psi_0=0.16\pi$$, $$k=3.525$$ Fig. 3 View largeDownload slide Magnitudes of coefficients for a skewed lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0.1,1.2]$$, $$a=0.005$$, and (a) $$\psi_0=0.15\pi$$, $$k=3.7$$, (b) $$\psi_0=0.16\pi$$, $$k=3.525$$ Figure 2(a) shows $$|A_j|$$ in a case where no Bloch waves are excited. These results were obtained by solving (60) by truncation at $$p=200$$ and $$q=200$$, and may be used as a benchmark for the effectiveness of the filtering method. The rate of decay in the coefficients is such that the error introduced by truncation is extremely small; in fact $$|A_j|$$ is numerically negligible for $$j\gtrapprox45$$. The conservation of energy criterion (105), and (102) with $$j=0$$ (the only propagating mode), are both satisfied to near machine precision. In Fig. 2(b), the frequency has been increased so that a single Bloch wave is excited, with $$\beta_y\approx1.78$$. No other parameters have been changed. Evidently, solving (60) directly by truncation would fail, since $$A_j\not\to0$$ as $$j\to\infty$$. However, the filtered coefficients $$A_j^{(1)}\to 0$$ extremely rapidly, and are numerically negligible for $$j\gtrapprox10$$. Computing $$A_j$$ again using infinite array subtraction produces results in agreement to $$13$$–$$14$$ significant figures. There is a slight drop in accuracy for the conservation of energy criterion, in that now only $$13$$ digits in $$E_1$$ cancel those in $$E_2$$. Similarly, (102) is now satisfied to $$14$$ significant figures. This is to be expected, in view of the more complex calculations that are needed when Bloch waves are excited. Figure 3 shows results for two more difficult cases. In Fig. 3(a), the lattice is skewed, and two Bloch waves are excited, with $$\beta_y\approx0.919$$, and $$\beta_y\approx1.67$$. Again, the filtered coefficients $$A_j^{(2)}$$ and the coefficients $$\hat{A}_j$$ decay very rapidly as $$j$$ is increased, but only to a magnitude around $$10^{-12}$$, so that the results are slightly less accurate in this case. This is confirmed by computing $$A_j$$ again, using infinite array subtraction; the results were found to be in agreement to $$12$$ significant figures. For these parameters, there are two propagating modes ($$j=0$$ and $$j=1$$), and (102) was satisfied to $$12$$ significant figures in both instances, with similar accuracy in the conservation of energy criterion (105). Finally, in Fig. 3(b), the same lattice is used, but the wavenumber and angle of incidence are chosen so that a single Bloch wave is excited with $$\beta_y\approx1.48$$, and a second is close to cutting on (there would be two Bloch waves in the equivalent plot for $$k=3.526$$). A Bloch wave close to cutting on corresponds to a value for $$\beta_y$$ with a small, positive imaginary part. In this case, terms up to $$j=500$$ and $$p=500$$ were retained in the calculations because the decay in $$A_j^{(1)}$$ and $$\hat{A}_j$$ is very slow. In fact this ceases after $${\sim}250$$ rows, when $$|A_j^{(1)}|\approx|\hat{A}_j|\approx10^{-12}$$. For some rows, the agreement between $$A_j$$ computed by filtering and the value obtained by infinite array subtraction is as low as $$10$$ significant figures, but (102) (with $$j=0$$ and $$j=1$$) and the conservation of energy criterion (105) are still satisfied to $$12$$ significant figures. 13. Concluding remarks We have demonstrated a new method for calculating the reflection and transmission properties of lattices, without the need for integral transforms or the Wiener–Hopf technique, and with no spurious reflection or damping effects. The method yields very accurate results, with slight deterioration close to band edges and when multiple Bloch waves are excited. The former issue could be countered by including Bloch vectors with small imaginary parts in the filtering process (see section 8). Several of the steps undertaken in setting up the new method are also needed if the same problem is to be solved via the Wiener–Hopf technique. In particular, the starting point is the reduction of a boundary value problem to an infinite linear system of equations using multipoles, and the same lattice sums and quasiperiodic Green’s functions arise automatically. Solving the propagation problem to determine the Bloch vectors for waves that can be excited inside the lattice is necessary for both approaches. For the Wiener–Hopf technique, these Bloch vectors correspond to zeroes of the kernel (or zeroes of the determinant of the kernel) that lie on the unit circle. Information about these is needed in order to construct the correct inversion contour (see (5) for details). The Wiener–Hopf method also requires factorisation of the kernel function into a product of functions analytic inside and outside the unit circle. For point scatterer problems in which a single amplitude coefficient is associated with each lattice element, this can be performed by numerically evaluating a Cauchy integral as in (7), or by approximating the field between consecutive rows as a finite expansion in grating modes and then locating kernel zeroes inside the unit circle as in (5). On the other hand, our method requires only the numerical solution of a rapidly convergent infinite linear system. The Wiener–Hopf method has one advantage, in that the analysis proves that the far field inside the lattice is composed of Bloch waves only. We have been unable to prove that this is the case using the current method, but we have verified it numerically. In fact the requirement that all other contributions must vanish provides a useful check on our results. The differences between the two approaches are much more stark when finite size effects are considered. In this case, each lattice element is associated with a sequence of amplitude coefficients, and the quasiperiodic Green’s functions and lattice sums required are those considered in (6, Appendices A and B). The Wiener–Hopf equation for such problems has a matrix kernel which cannot be factorised by any known method. The approximate scheme used from (6) can be used instead, but this is complicated and very difficult to implement. On the other hand, our method generalises immediately. The algebra becomes more involved, but the procedure for obtaining solutions is largely identical to the point scatterer case. Semi-infinite lattices composed of finite-sized scatterers will be considered in the second part of this work. The technique can also be generalised to account for aperiodic sources (such as line sources) and lattices with defects, by using the array scanning method (22, 23), an approach which was developed for linear arrays in (24) and (25). References 1. Nicorovici N. A. McPhedran R. C. and Botten L. C. Photonic band gaps for arrays of perfectly conducting cylinders, Phys. Rev. E  52 ( 1995) 1135– 1145. Google Scholar CrossRef Search ADS   2. Sigalas M. and Economou E. N. Elastic and acoustic-wave band-structure, J. Sound Vib.  158 ( 1992) 377– 382. Google Scholar CrossRef Search ADS   3. Sigalas M. and Economou E. N. Band structure of elastic waves in two-dimensional systems, Solid State Commun.  86 ( 1993) 141– 143. Google Scholar CrossRef Search ADS   4. Joannopoulos J. D. Johnson S. G. Winn J. N. and Meade R. D. Photonic Crystals. Molding the Flow of Light , 2nd edn. ( Princeton University Press, Princeton, USA and Oxford, UK 2008). 5. Tymis N. and Thompson I. Low-frequency scattering by a semi-infinite lattice of cylinders, Q. J. Mech. Appl. Math.  64 ( 2011) 171– 195. Google Scholar CrossRef Search ADS   6. Tymis N. and Thompson I. Scattering by a semi-infinite lattice of cylinders and the excitation of Bloch waves, Q. J. Mech. Appl. Math.  67 ( 2014) 469– 503. Google Scholar CrossRef Search ADS   7. Haslinger S. G. Craster R. V. Movchan A. B. Movchan N. V. and Jones. I. S. Dynamic interfacial trapping of flexural waves in structured plates. Proc. R. Soc. A  472 ( 2016) 20150658 ( 25 pages). Google Scholar CrossRef Search ADS   8. Martin P. A. Multiple Scattering. Interaction of Time-Harmonic Waves with N Obstacles . ( Cambridge University Press, Cambridge, UK 2006). Google Scholar CrossRef Search ADS   9. Hills N. L. and Karp S. N. Semi-infinite diffraction gratings—I, Comm. Pure Appl. Maths  18 ( 1965) 203– 233. Google Scholar CrossRef Search ADS   10. Hills N. L. Semi-infinite diffraction gratings. II. Inward resonance, Comm. Pure Appl. Maths  18 ( 1965) 389– 395. Google Scholar CrossRef Search ADS   11. Linton C. M. and Martin P. A. Semi-infinite arrays of isotropic point scatterers. A unified approach, SIAM J. Appl. Math.  64 ( 2004) 1035– 1056. Google Scholar CrossRef Search ADS   12. Noble B. Methods Based on the Wiener–Hopf Technique . ( Chelsea, New York 1988). 13. Lawrie J. B. and Abrahams I. D. A brief historical perspective of the Wiener–Hopf technique. J. Engng. Math.  59 ( 2007) 351– 358. Google Scholar CrossRef Search ADS   14. Rogosin S. and Mishuris G. Constructive methods for factorization of matrix-functions, IMA J. Appl. Math.  81 ( 2016) 365– 391. Google Scholar CrossRef Search ADS   15. Linton C. M. Porter R. and Thompson I. Scattering by a semi-infinite periodic array and the excitation of surface waves, SIAM J. Appl. Math.  67, ( 2007) 1233– 1258. Google Scholar CrossRef Search ADS   16. Foldy L. L. The multiple scattering of waves I. General theory of isotropic scattering by randomly distributed scatterers, Phys. Rev.  67 ( 1945) 107– 119. Google Scholar CrossRef Search ADS   17. Brillouin L. Wave Propagation in Periodic Structures , 2nd edn. ( Dover, New York 1953). 18. Linton C. M. and Thompson I. Resonant effects in scattering by periodic arrays, Wave Motion  44 ( 2007) 167– 175. Google Scholar CrossRef Search ADS   19. Botten L. C. Nicorovici N. A. Asatryan A. A. McPhedran R. C. Martijn de Sterke C. and Robinson P. A. Formulation for electromagnetic scattering and propagation through grating stacks of metallic and dielectric cylinders for photonic crystal calculations. Part II. Properties and implementation. J. Opt. Soc. Am. A  17 ( 2000) 2177– 2190. Google Scholar CrossRef Search ADS   20. Twersky V. On the scattering of waves by an infinite grating, IRE Trans. on Antennas and Propagation  4 ( 1956) 330– 345. Google Scholar CrossRef Search ADS   21. Linton C. M. Lattice sums for the Helmholtz equation, SIAM Review  52 ( 2010) 630– 674. Google Scholar CrossRef Search ADS   22. Wu C. P. and Galindo V. Properties of a phased array of rectangular waveguides with thin walls, IEEE Trans. Antennas Propagat.  14 ( 1966) 163– 173. Google Scholar CrossRef Search ADS   23. Munk B. A. and Burrell G. A. Plane-wave expansion for arrays of arbitrarily oriented piecewise linear elements and its application in determining the impedance of a single linear antenna in a lossy half-space, IEEE Trans. Antennas Propagat.  27 ( 1979) 331– 343. Google Scholar CrossRef Search ADS   24. Thompson I. and Linton C. M. On the excitation of a closely spaced array by a line source, IMA J. Appl. Math.  72 ( 2007) 476– 497. Google Scholar CrossRef Search ADS   25. Thompson I. and Linton C. M. An interaction theory for scattering by defects in arrays, SIAM J. Appl. Math.  68 ( 2008) 1783– 1806. Google Scholar CrossRef Search ADS   Appendix Consider the function $$f$$, defined via   $$\label{def:f} f(z_1,\ldots,z_{\lambda+1}) = \sum_{m=1}^\lambda \prod_{\substack{n=1 \\ n \neq m}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_m}.$$ (A.1) We will prove that this is in fact the constant function $$f(z_1,\ldots,z_{\lambda+1})\equiv 1$$. First consider $$f$$ as a function of the complex variable $$z_1$$, with $$z_2,\ldots,z_{\lambda+1}$$ fixed (and distinct). For any integer $$p$$ between $$2$$ and $$\lambda$$ (inclusive), there are two ways in which $$z_1-z_p$$ can appear in the denominator: either $$m=1$$ or $$m=p$$. Separating these terms from the sum, and considering the limit in which $$z_1\to z_p$$, we find that   \begin{align*} \lim_{z_1\to z_p} \Biggl[ \prod_{n=2}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_1} + \prod_{\substack{n=1 \\ n \neq p}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_p} \Biggr] &= \prod_{\substack{n=2 \\ n \neq p}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_p} \lim_{z_1\to z_p} \Biggl[ \frac{z_p-z_{\lambda+1}}{z_p-z_1} + \frac{z_1-z_{\lambda+1}}{z_1-z_p} \Biggr] \\ &= \prod_{\substack{n=2 \\ n \neq p}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_p}. \end{align*} Therefore, $$f$$ is an entire function of the complex variable $$z_1$$. Furthermore, we may observe that   $$f(z_1,\ldots,z_{\lambda+1}) = \prod_{n=2}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_1} + \sum_{m=2}^\lambda \frac{z_1-z_{\lambda+1}}{z_1-z_m} \prod_{\substack{n=2 \\ n \neq m}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_m},$$ (A.2) from which we immediately find that   $$\lim_{z_1\to\infty } f(z_1,\ldots,z_{\lambda+1}) = \sum_{m=2}^\lambda \prod_{\substack{n=2 \\ n \neq m}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_m}.$$ (A.3) Hence, $$f$$ is a bounded entire function of $$z_1$$, so Liouville’s theorem shows that it must be independent of $$z_1$$. The same calculation can be made with $$z_n$$ in place of $$z_1$$, for $$n=2,\ldots,\lambda$$. Since replacing $$z_n$$ with $$z_n-\alpha$$ for $$n=1,\ldots,\lambda$$ does not affect $$f$$, we may make these translations in (A.1) to show that it is also independent of $$z_{\lambda+1}$$. Finally, we set $$z_1=z_{\lambda+1}$$, so that the terms with $$m>1$$ disappear from (A.1), and we easily find that   $$f(z_{\lambda+1},z_2,\ldots,z_{\lambda+1}) = 1,$$ (A.4) which completes the proof. © The Author, 2017. Published by Oxford University Press; all rights reserved. For Permissions, please email: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png The Quarterly Journal of Mechanics and Applied Mathematics Oxford University Press

# A Direct Method for Bloch Wave Excitation by Scattering at the Edge of a Lattice. Part I: Point Scatterer Problem

, Volume 71 (1) – Feb 1, 2018
24 pages

/lp/ou_press/a-direct-method-for-bloch-wave-excitation-by-scattering-at-the-edge-of-w7lsHwMOdl
Publisher
Oxford University Press
ISSN
0033-5614
eISSN
1464-3855
D.O.I.
10.1093/qjmam/hbx022
Publisher site
See Article on Publisher Site

### Abstract

Summary A new method for determining the reflection and transmission properties of lattices is developed. The method uses multipole expansions, and certain transformations of the algebraic equation systems that appear when boundary conditions are applied. It is more direct, and much simpler, than earlier approaches based on integral transforms and the Wiener–Hopf technique. The method is demonstrated for the case of a semi-infinite lattice of sound soft acoustic point scatterers, but can easily be generalised to account for finite size effects, and more general boundary conditions. 1. Introduction Bloch waves propagate through periodic structures without loss of energy, and have been found to exist in a number of physical contexts including electromagnetism, acoustics and elasticity (1, 2, 3). Generally, a periodic structure that supports Bloch wave propagation can only do so at certain frequencies, and in certain directions. This band structure gives rise to a range of applications, including basic components such as filters and mirrors, and more radical ideas such as optical circuits (4). Much of the literature in this area is devoted to the propagation problem, that is determining the parameter regimes in which Bloch waves can exist in a particular medium. In contrast, very little attention has been paid to the excitation problem, in which a wave incident from free space strikes the edge of a periodic structure, and some of the energy may be converted into Bloch waves. The excitation problem is generally much more difficult than the associated propagation problem, but it has enormous implications for the effectiveness of periodic media in manipulating wave propagation. For example, a band gap filter can prevent the transmission of signals at certain frequencies, but part of the signal at the desired frequencies may also be reflected back. By solving the excitation problem, it is possible to determine the proportion of incident energy that is lost in this way. A qualitative discussion of electromagnetic Bloch wave excitation is given in (4, pp. 221–227), but here the authors are only concerned with the number of diffracted Bloch waves and the directions in which these propagate. No method for determining their amplitudes is given. Perhaps the simplest case in which Bloch wave excitation can occur involves a plane wave incident on a semi-infinite lattice of acoustic point scatterers. This was considered in (5). In a subsequent paper, the same authors went on to introduce finite size effects, thereby allowing for a greater range of boundary conditions and frequencies (6). Another recent paper considers excitation of Bloch waves in a semi-infinite lattice of pinned points on a thin plate modelled by Kirchhoff theory (7). In all these three papers, the boundary value problem is converted into an infinite linear system of algebraic equations (using the multipole method (8, Chapter 4), or a reduced form for point scatterers), but these systems can only be solved directly in cases where no Bloch waves exist; otherwise truncation introduces spurious reflection effects which change the nature of the problem. However, applying an appropriate integral transform leads to a Wiener–Hopf equation. Once this is solved, the reflection and transmission properties of the lattice can be calculated. The idea of using the Wiener–Hopf method in this way was originally introduced in a study of wave diffraction by a linear array (9, 10) (see also (11)). For problems involving point scatterers, the Wiener–Hopf equation is scalar, and can be solved by a standard procedure (12). On the other hand, allowing for finite size effects leads to a matrix Wiener–Hopf equation. Although solutions to equations of this type are known to exist, obtaining them is notoriously difficult, and a general method has yet to be found (see (13) for a discussion of this issue and relevant references). In particular, the matrix equation obtained in (6) is not amenable to any of the methods recently surveyed in (14). The authors of (6) overcame this issue by approximating the field between consecutive rows as a finite expansion in grating modes. This eventually reduces the problem of solving the Wiener–Hopf equation to a complex root finding problem. However, implementing a solution to this turns out to be rather difficult, because the roots must be located with great precision, and the complexity of the function and the number of roots become greater as the frequency is increased. The purpose of this article is to introduce a more direct method for calculating the reflection and transmission properties of lattices, based on the filtering technique introduced in (15). The idea behind this is very straightforward. First, the dispersion relation from the associated propagation problem is solved, to determine the Bloch vectors for waves that may be excited inside the lattice. The algebraic equations obtained from the multipole method are then transformed in such a way that slowly convergent terms, present due to the excitation of Bloch waves, are removed. Only Bloch vectors corresponding to waves that transport energy towards infinity are included in this process; this is the radiation condition for Bloch wave excitation problems. The resulting rapidly convergent system can be truncated without introducing spurious reflection or damping effects. After inverting the truncated system numerically, the solution to the original system can be reconstructed, and all important quantities such as the Bloch wave amplitudes, and the proportions of energy transmitted into and reflected back from the lattice, can be calculated. Aside from being simpler than the Wiener–Hopf approach discussed above, our new method has the advantage that allowing for lattice elements of finite size has very little effect aside from introducing more complicated algebra. Therefore it makes sense to present the method initially for the case of acoustic point scatterers. The structure of the article is as follows. Sections 2–4 contain a review of the basic material needed for modelling wave interactions with lattices. Section 5 contains some formulas for the evaluation of quasiperiodic Green’s functions, which are required by our method. The same Green’s functions are used when the problem is solved via the Wiener–Hopf technique (see (5)), so the need for these is not a new complication. Section 6 presents a simple method for solving the propagation problem, using the quasiperiodic Green’s functions from section 5. The algebraic system of equations that governs scattering by a semi-infinite lattice is obtained in section 7. This is the same equation that was solved via the Wiener–Hopf method in (5). Indeed, all of the material in sections 2–7 can be found in the references given. Its inclusion here is necessary in order to present our method in a coherent fashion, using consistent notation. The application of the filtering method to the infinite algebraic system is performed in section 8. This is the key part of the article, and the process for transforming the system in the presence of an arbitrary number of Bloch waves is the main result. To maintain the flow of the material, a technical part of the proof that establishes the validity of the transformation is deferred to the appendix. Section 9 contains details of another idea from (15), namely infinite array subtraction. This can be used to recalculate certain coefficients in the problem, given the amplitudes of the excited Bloch waves obtained from the filtering method, and provides a stringent means of testing the accuracy of our results. Formulae for determining the reflected and transmitted fields, and a related conservation of energy criterion are obtained in sections 10 and 11. A small number of numerical results are presented in section 12. The purpose of these is to demonstrate that the method works, and to illustrate its accuracy. A greater range of results will be presented in the second part of this work, where finite size effects will be introduced. Some concluding remarks, including a detailed comparison of the filtering method and the Wiener–Hopf approach, are given in section 13. 2. Scattering by a small cylinders We will consider the interaction of two-dimensional wave fields $$u(\mathbf{r})$$ with cylinders of radius $$a$$ and infinite length (so that there is no variation in the $$z$$ direction). Outside the cylinders, the field $$u(\mathbf{r})$$ must satisfy the Helmholtz equation   $$\label{helmholtz} (\nabla^2+k^2)u(\mathbf{r}) = 0,$$ (1) and on their surfaces we have the Dirichlet boundary condition   $$\label{dirichlet} u(\mathbf{r}) = 0.$$ (2) After $$u(\mathbf{r})$$ has been determined, a time-harmonic acoustic potential can be retrieved by writing   $$U(\mathbf{r},t) = \mathrm{Re}[ u(\mathbf{r}){\mathrm e}^{-{\mathrm i}\omega t}],$$ (3) where $$\omega=ck$$, with $$\omega$$ representing frequency, and $$c$$ the speed of sound in the exterior medium. Note that $$U(\mathbf{r},t)$$ can also be interpreted as the $$z$$ component of a transverse magnetic (TM) electric field vector incident on a lattice of perfectly conducting cylinders. In this case $$c$$ represents the speed of light. Now any field incident on a cylinder centred at $$\mathbf{r}=\mathbf{R}$$ can be expanded in the form   $$u^{\mathrm i}(\mathbf{r}) = \sum_{n=-\infty}^\infty I_n {J}_n(k|\mathbf{r}-\mathbf{R}|){\mathrm e}^{{\mathrm i} n\mu},$$ (4) where $${J}_n$$ is the Bessel function of order $$n$$, and $$\mu$$ is the anticlockwise angle between the positive $$x$$ axis and the vector $$\mathbf{r}-\mathbf{R}$$. Setting $$\mathbf{r}=\mathbf{R}$$ reveals that $$I_0=u^{\mathrm i}(\mathbf{R})$$, and hence, on the surface of the cylinder we have   $$\label{expand_inc} u^{\mathrm i}(\mathbf{r}) = u^{\mathrm i}(\mathbf{R}) {J}_0(ka) + O(ka).$$ (5) If we take the scattered field to be   $$\label{single_scattered} u^{{\mathrm s}}(\mathbf{r}) = A {\mathrm{H}}_0^{(1)}(k|\mathbf{r}-\mathbf{R}|),$$ (6) where $${\mathrm{H}}_0^{(1)}$$ is a Hankel function of the first kind, then the Dirichlet boundary condition (2) is satisfied at the leading-order if the amplitude coefficient $$A$$ is given by   $$\label{foldy} A = -Z_0 u^{\mathrm i}(\mathbf{R}),$$ (7) in which   $$\label{def:Z0} Z_0 = \frac{{J}_0(ka)}{{\mathrm{H}}_0^{(1)}(ka)} = \biggl( 1 + \frac{2{\mathrm i}}{\pi}\bigl[ C + \ln(ka) - \ln 2\bigr]\biggr)^{-1} + O\bigl((ka)^2\bigr).$$ (8) Here, $$C=0.5772\ldots$$ is Euler’s constant. This method for treating scattering in the low-frequency limit dates back to (16). The same results can also be obtained using asymptotic expansions (see (8, Chapter 8)). In this way one can also derive approximate scattering coefficients for arbitrary small shapes. Higher-order terms can be obtained using shift operators (8, Sections 2.4–2.5), leading to Graf’s addition theorem and the full linear theory for multiple scattering. However, the leading-order approximation, in which $$O(ka)$$ and smaller terms are discarded, is sufficient for our purposes. 3. Two-dimensional lattices A two-dimensional lattice can be formed by placing identical objects at points with position vectors   $$\label{def:Rjp} \mathbf{R}_{jp} = j\mathbf{s}_1 + p\mathbf{s}_2, \quad j,p\in\mathbb{Z}.$$ (9) Without loss of generality, we may assume that   $$\label{def:s12} \mathbf{s}_1 = s_1 \hat{\mathbf{x}} \quad\text{and}\quad \mathbf{s}_2 = \eta_1 \hat{\mathbf{x}} + \eta_2 \hat{\mathbf{y}},$$ (10) where the circumflex denotes a unit vector,   $$s_1/2 \ge \eta_1 \ge 0 \quad\text{and}\quad \eta_2>0.$$ (11) Here we have introduced the convention that $$|\mathbf{v}|=v$$ for any vector $$\mathbf{v}$$, which will be used throughout. We will also require the reciprocal lattice vectors   $$\label{def:recip_vecs} \mathbf{s}_1^* = \frac{1}{s_1} \left(\hat{\mathbf{x}}-\frac{\eta_1}{\eta_2}\,\hat{\mathbf{y}}\right) \quad\text{and}\quad \mathbf{s}_2^* = \frac{\hat{\mathbf{y}}}{\eta_2}.$$ (12) These have the property that   $$\label{recip_prop} \mathbf{s}_j \cdot \mathbf{s}_p^* = \delta_{jp},$$ (13) where $$\delta_{jp}$$ is the Kronecker delta. Consequently, if   $$\label{def:Rmn} \mathbf{R}_{mn}^* = 2\pi(m \mathbf{s}_1^* + n \mathbf{s}_2^*), \quad m,n\in\mathbb{Z},$$ (14) then   $$\label{recip_prop2} {\mathrm e}^{{\mathrm i} \mathbf{R}_{jp}\cdot \mathbf{R}_{mn}^*} = 1,$$ (15) for all combinations of integers $$j$$, $$p$$, $$m$$ and $$n$$. Further details of lattice theory can be found in (17) and (4, Appendix B). 4. Grating modes In what follows, we will consider wave fields that have the one-dimensional quasiperiodicity property   $$\label{qp_1d} u(\mathbf{r} + s_1\hat{\mathbf{x}} ) = {\mathrm e}^{{\mathrm i} s_1 \beta_x} u(\mathbf{r}).$$ (16) A field that satisfies (16) and the Helmholtz equation (1) may be expanded as a series of grating modes. That is,   $$\label{grating_expansion} u(\mathbf{r}) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} \beta_{xj} x} \bigl[ c_j^- {\mathrm e}^{\gamma(\beta_{xj})y} + c_j^+ {\mathrm e}^{-\gamma(\beta_{xj})y} \bigr],$$ (17) where   $$\label{def:bxj} \beta_{xj} =\beta_x+2j\pi/s_1, \quad j\in\mathbb{Z},$$ (18) and   $$\label{def:gamma} \gamma(t) = \begin{cases} \sqrt{ t^2 - k^2 } &\text{if } |t|\ge k, \\ -{\mathrm i} \sqrt{ k^2 - t^2 } &\text{if } |t|< k. \end{cases}$$ (19) In general, the amplitude coefficients $$c_j^\pm$$ are different in different regions, for example, above and below a grating, or between different rows of a lattice. There are a finite number terms in (17) for which the real part of $$\gamma(\beta_{xj})$$ is zero; these represent propagating modes. Since the distinction between propagating modes and evanescent modes is important, we introduce the sets   $$\label{def:MN} \mathcal{M} = \bigl\{ j: \ j\in\mathbb{Z}, \ |\beta_{xj}|\le k \} \quad\text{and}\quad \mathcal{N} = \bigl\{ j: \ j\in\mathbb{Z}, \ |\beta_{xj}| > k \}.$$ (20) We will refer to modes with amplitude coefficient $$c_j^+$$ as upward oriented (propagating in the positive $$y$$ direction, or decaying as $$y\to\infty$$), whereas those with coefficients $$c_j^-$$ are downward oriented. To avoid certain technicalities associated with Wood’s anomalies (also sometimes called resonances), we will assume that the parameters $$\beta_x$$, $$s_1$$ and $$k$$ are such that $$|\beta_{xj}|\neq k$$ for any integer $$j$$, so that $$\gamma(\beta_{xj})\neq0$$. For the case of a linear array, a method for dealing with Wood’s anomalies was developed in (18). It is possible to adapt the analysis of semi-infinite lattices in a similar way, but this is not pertinent to the current article. For later purposes, it is useful to determine the direction in which a field of the form (17) transports energy across lines parallel to the $$x$$ axis. Now the average energy flux over any line $$\mathcal S$$ in one time-period is given by the line integral (5)   $$\label{energy_int} \langle E_{\mathcal S} \rangle = -\,\frac{P_0\omega}{2} \mathrm{Im} \int_{\mathcal S} u(\mathbf{r}) \frac{ \partial }{\partial n} u^*(\mathbf{r}) \, {\mathrm d} s,$$ (21) where the derivative is taken in a direction normal to $$\mathcal S$$ and the superscript ‘$$*$$’ denotes a complex conjugate. If $$\langle E_{\mathcal S} \rangle > 0$$, the direction of net energy flux is the same as the orientation of the normal. The constant $$P_0$$ represents the quiescent fluid density in the case of an acoustic wavefield. For a TM polarised electric field, this should be replaced by $$1/(\omega^2\mu_0)$$, where $$\mu_0$$ is the magnetic permeability of the exterior medium. Taking $$\mathcal S$$ to be the straight line from the point $$(x_0 - s_1/2,y_0)$$ to $$(x_0 + s_1/2,y_0)$$, and directing the normal parallel to $$\hat{\mathbf{y}}$$, (21) becomes   $$\label{energy2} \langle E_{\mathcal S} \rangle = -\,\frac{P_0\omega}{2} \mathrm{Im} \int_{x_0-s_1/2}^{x_0+s_1/2} u(\mathbf{r}) \frac{ \partial }{\partial y} u^*(\mathbf{r})\biggr|_{y=y_0} \, {\mathrm d} x.$$ (22) Now   $$\label{grating_int} \int_{x_0-s_1/2}^{x_0+s_1/2} {\mathrm e}^{{\mathrm i} \beta_{xj}x} {\mathrm e}^{-{\mathrm i} \beta_{xp}x} \, {\mathrm d} x = s_1 \delta_{jp},$$ (23) so if we use (17) in (22), we obtain   \begin{multline} \langle E_{\mathcal S} \rangle = -\,\frac{P_0\omega s_1}{2} \mathrm{Im} \sum_{j=-\infty}^\infty \gamma^*(\beta_{xj}) \bigl[ c_j^- {\mathrm e}^{\gamma(\beta_{xj})y_0} + c_j^+ {\mathrm e}^{-\gamma(\beta_{xj})y_0} \bigr] \\ \times \bigl[ (c_j^-)^* {\mathrm e}^{\gamma^*(\beta_{xj})y_0} - (c_j^+)^* {\mathrm e}^{-\gamma(\beta_{xj})^*y_0} \bigr]. \end{multline} (24) Considerable simplifications now occur on separating terms in which $$\gamma(\beta_{xj})$$ is imaginary from those in which it is real. We find that   $$\label{energy_prop} \langle E_{\mathcal S} \rangle = \frac{P_0\omega s_1}{2} \sum_{j\in\mathcal M} |\gamma(\beta_{xj})| \bigl( |c_j^+|^2 - |c_j^-|^2 \bigr) - P_0\omega s_1 \sum_{j\in\mathcal N}\gamma(\beta_{xj})\mathrm{Im}\bigl[ c_j^+ (c_j^-)^*\bigr],$$ (25) where the sets $$\mathcal{M}$$ and $$\mathcal{N}$$ are defined in (20). The same formula is given in (6) and also earlier (with different notation) in (19). As noted in (19), the second term on the right-hand side of (25) plays no role in problems involving a single grating illuminated by a plane wave, because in such cases there are no incoming evanescent modes, so for $$j\in\mathcal N$$ we have $$c_j^-=0$$ above the grating and $$c_j^+=0$$ below. However, this term is of crucial importance in lattice problems because there is a full set of upward- and downward-oriented grating modes between each consecutive pair of rows. Since $$x_0$$ and $$y_0$$ do not appear in (25), there is some freedom in setting the location of $$\mathcal S$$. Convenient choices are typically $$x_0=q\eta_1$$ and $$y_0=(q+1/2)\eta_2$$ for an integer $$q$$, so that there are no singularities on the integration path (see, for example, Fig. 7 in (5)). 5. Quasiperiodic Green’s functions The quasiperiodic Green’s function for a row of sources located at the points $$\mathbf{r}=\mathbf{R}_{j0}$$ (9) is given by   $$\label{def:QPG1d} G_0(\mathbf{r};\beta_x) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} js_1 \beta_x} {\mathrm{H}}_0^{(1)}(k|\mathbf{r}-\mathbf{R}_{j0}|).$$ (26) By a straightforward application of Poisson summation (20, 5), this may be rewritten in the spectral form   $$\label{spectral:QPG1d} G_0(\mathbf{r};\beta_x) = -\,\frac{2{\mathrm i}}{s_1}\sum_{j=-\infty}^\infty \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x} }{\gamma(\beta_{xj})}\, {\mathrm e}^{- \gamma(\beta_{xj}) |y| },$$ (27) where $$\beta_{xj}$$ and $$\gamma$$ are defined in (18) and (19), respectively. Closely related to $$G_0$$ is the one-dimensional lattice sum   \begin{align} \sigma_0(\beta_x) &= \lim_{r\to 0} \bigl[ G_0(\mathbf{r};\beta_x) - {\mathrm{H}}_0^{(1)}(kr) \bigr] \\ \end{align} (28)  \begin{align} &= 2\sum_{j=1}^\infty \cos(js_1\beta_x){\mathrm{H}}_0^{(1)}(js_1), \label{def:sigma0} \end{align} (29) which is sometimes called a Schlömilch series, and can be evaluated using the methods in (21). We will also need the quasiperiodic Green’s functions for multiple rows, which were originally considered in (5). For rows $$q_0,\ldots,q_1$$, we have   \begin{align} G_0^{(q_0,q_1)}(\mathbf{r};{\boldsymbol{\beta}}) &= \sum_{q=q_0}^{q_1} \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} \mathbf{R}_{jq} \cdot {\boldsymbol{\beta}} } {\mathrm{H}}_0^{(1)}(k|\mathbf{r}-\mathbf{R}_{jq}|) \label{def:22QPG} \\ \end{align} (30)  \begin{align} &= \sum_{q=q_0}^{q_1} {\mathrm e}^{{\mathrm i} q \mathbf{s}_2\cdot{\boldsymbol{\beta}}} G_0(\mathbf{r}-q\mathbf{s}_2;\beta_x), \label{22QPG2} \end{align} (31) where $${\boldsymbol{\beta}}=\beta_x\hat{\mathbf{x}}+\beta_y\hat{\mathbf{y}}$$. If $$q_0$$ and $$q_1$$ are finite integers, this has the spectral form   \begin{align} G_0^{(q_0,q_1)}(\mathbf{r};{\boldsymbol{\beta}})= -\,\frac{2{\mathrm i}}{s_1}\sum_{j=-\infty}^\infty \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x \mp \gamma(\beta_{xj})y } }{\gamma(\beta_{xj})}\,\frac{ {\mathrm e}^{q_0 w_j^{^\pm}} - {\mathrm e}^{(1+q_1)w_j^{^\pm}} } { 1 - {\mathrm e}^{w_j^{^\pm}} }, \label{MRQPG} \end{align} (32) where   \begin{align} w^{\pm}_j = \pm\eta_2 \gamma(\beta_{xj}) + {\mathrm i}( \eta_2\beta_y - 2j\pi\eta_1/s_1), \label{def:wj} \end{align} (33) and the upper and lower signs are to be taken when $$y\ge q_1\eta_2$$ and $$y\le q_0\eta_2$$, respectively. The situation in which $$q_0\eta_2 < y < q_1\eta_2$$ can be accounted for by considering the rows above and below the observer separately and using (32) twice. Finally, for semi-infinite lattices, (5) gives the formulae   \begin{align} G_0^{(q_0,\infty)}(\mathbf{r};{\boldsymbol{\beta}}) = \frac{2{\mathrm i}}{s_1} \sum_{j=-\infty}^\infty {\mathrm e}^{q_0 w_j^{^-}} \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x + \gamma(\beta_{xj})y } }{\gamma(\beta_{xj}) \big({\mathrm e}^{w_j^{^-}}-1 \big) }, \quad y \le q_0\eta_2 \label{QPG:upper} \end{align} (34) and   \begin{align} G_0^{(-\infty,q_1)}(\mathbf{r};{\boldsymbol{\beta}}) = \frac{2{\mathrm i}}{s_1} \sum_{j=-\infty}^\infty {\mathrm e}^{q_1 w_j^{^+}} \frac{{\mathrm e}^{{\mathrm i} \beta_{xj} x -\gamma(\beta_{xj})y } }{\gamma(\beta_{xj}) \big({\mathrm e}^{-w_j^{^+}}-1\big) }, \quad y \ge q_1\eta_2. \label{QPG:lower} \end{align} (35) 6. The propagation problem The problem of Bloch wave propagation through two-dimensional lattices of cylindrical Dirichlet scatterers was considered in detail in (1). Here we require only a leading-order approximation, valid when the scatterer radius is much smaller than the incident wavelength. In this case a Bloch wave has the form   $$\label{Bloch1} u({\mathbf{{r}}}) = B G_0^{(-\infty,\infty)}({\mathbf{{r}}};{\boldsymbol{\beta}}),$$ (36) for a Bloch vector $${\boldsymbol{\beta}}$$ and an arbitrary amplitude coefficient $$B$$. This has the two-dimensional quasiperiodicity property   $$\label{QP2d} u({\mathbf{{r}}}+{\mathbf{{R}}}_{jp}) = {\mathrm e}^{{\mathrm i} {\mathbf{{R}}}_{jp}\cdot {\boldsymbol{\beta}}} u({\mathbf{{r}}}),$$ (37) which incorporates the one-dimensional property (16). The Bloch vector must be chosen so that the boundary condition on the cylinder centred at the origin is satisfied. Boundary conditions elsewhere then follow by quasiperiodicity. From the perspective of this cylinder, the Bloch wave may be decomposed into incoming and scattered fields $$u_0^{\mathrm i}$$ and $$u_0^{\mathrm s}$$, where the incoming field consists of the radiation from all the other cylinders. That is,   \begin{align} u^{\mathrm i}_0({\mathbf{{r}}};{\boldsymbol{\beta}}) = B\sum_{p=-\infty}^\infty \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} {\mathbf{{R}}}_{jp}\cdot {\boldsymbol{\beta}}} {H}_0^{(1)}(k|{\mathbf{{r}}}- {\mathbf{{R}}}_{jp}|) \label{def:u0i} \end{align} (38) and   \begin{align} u^{\mathrm s}_0({\mathbf{{r}}};{\boldsymbol{\beta}}) = B {H}_0^{(1)}(kr), \end{align} (39) where the prime in (38) indicates that the term in which $$R_{jp}=0$$ is to be omitted from the summation. Comparing these to (6) and (7), we find that $$B=-Z_0 u_0^{\mathrm i}({\mathbf{{0}}};{\boldsymbol{\beta}})$$, where $$Z_0$$ is given by (8). Evaluating (38) at the origin then yields   $$\label{bloch_dispersion} 1 + Z_0 \Xi_0({\boldsymbol{\beta}}) = 0,$$ (40) where   $$\label{def:2dsigma} \Xi_0({\boldsymbol{\beta}}) = \sum_{p=-\infty}^\infty \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} {\mathbf{{R}}}_{jp}\cdot {\boldsymbol{\beta}}} {H}_0^{(1)}(k R_{jp}).$$ (41) A number of methods for computing two-dimensional lattice sums of this type are discussed in (21). One simple approach is to use the decomposition   $$\label{2dsigma_decomp} \Xi_0({\boldsymbol{\beta}}) = G_0^{(-\infty,-1)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + \sigma_0(\beta_x) + G_0^{(1,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}),$$ (42) where the terms on the right-hand side are defined in section 5. Equation (40) is the dispersion relation for Bloch waves. Given fixed values for all parameters except $$\beta_y$$, its real roots must be determined before we can proceed to consider Bloch wave excitation. To simplify this process, we use results from (21, Section 3) to conclude that $${\text{Re}}[\Xi_0]=-1$$, and so (40) becomes   $$\label{d2} 1 + W_0\,{\text{Im}}[\Xi_0({\boldsymbol{\beta}})] = 0,$$ (43) where   $$W_0 = {\mathrm i} Z_0/(1-Z_0) = {J}_0(ka) / {Y}_0(ka).$$ (44) Here, $${Y}_0$$ is the Neumann function of order zero, and the coefficient $$W_0$$ is clearly real. Now, from (9) and (10), we have $${\mathbf{{R}}}_{jp}\cdot\hat{{\mathbf{{y}}}}=p\eta_2$$, so (37) shows that all possibilities are accounted for if we consider   $$\label{beta_y_range} 0 \le \beta_y <2\pi/\eta_2.$$ (45) We must also take into account the fact that the function $$\Xi_0$$ has poles at points where $$|{\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*| = k$$, for some $$m,n\in\mathbb{Z}$$ (21, (3.18)). From (14) and (12),   $$\hat{{\mathbf{{x}}}} \cdot ({\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*) = \beta_x + \frac{2m\pi}{s_1}$$ (46) so at any pole we must have   $$\label{m_range} \left|\beta_x + \frac{2m\pi}{s_1} \right| \le k.$$ (47) Given a value for $$m$$ which satisfies (47), there are two possible values for the $$y$$ component. These are given by   $$\hat{{\mathbf{{y}}}}\cdot ({\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*) = \pm\sqrt{k^2- \bigl(\hat{{\mathbf{{x}}}}\cdot ({\boldsymbol{\beta}} + {\mathbf{{R}}}_{mn}^*)\bigr)^2 },$$ (48) from which we easily obtain   $$\beta_y + \frac{2n\pi}{\eta_2} = \frac{2\pi m\eta_1}{s_1 \eta_2} \pm \sqrt{k^2 - (\beta_x+2m\pi/s_1)^2 }.$$ (49) Since $$n\in\mathbb{Z}$$ it now follows that, for each $$m$$ satisfying (47), there are precisely two poles such that $$\beta_y$$ satisfies (45), unless the square root evaluates to zero, in which case there is only one pole for this particular choice of $$m$$. Another useful piece of information that can be obtained from the propagation problem is the direction in which each Bloch wave carries energy across lines parallel to the $$x$$ axis. Now Bloch waves do not have well-defined phase velocities, because integer multiples of $${\mathbf{{s}}}_1^*$$ and $${\mathbf{{s}}}_2^*$$ can be added to $${\boldsymbol{\beta}}$$ with no substantive effect in view of (15) and (37). Therefore, we must use information about energy transportation to determine the direction of propagation (see (4, pp. 222–225) for further discussion). To this end, we evaluate (25) for the Bloch wave represented by (36). Thus, using (34) and (35) with $$q_1=q$$, and $$q_0=1+q$$ we find that for $$q \eta_2 \le y \le (1+q) \eta_2$$, a Bloch wave can be expanded in the form (17), with   $$c_j^+ = \frac{2{\mathrm i} B {\mathrm e}^{q w_j^+}}{s_1\gamma(\beta_{xj})( {\mathrm e}^{-w_j^+} - 1 )} \quad\text{and}\quad c_j^- = \frac{2{\mathrm i} B {\mathrm e}^{(1+q) w_j^-}}{s_1\gamma(\beta_{xj})( {\mathrm e}^{w_j^-} - 1 )}.$$ (50) If $$\gamma(\beta_{xj})$$ is imaginary, then so too are $$w_j^+$$ and $$w_j^-$$, and so   $$\label{cjsq} |c_j^+|^2 = \frac{4 |B|^2}{s_1^2\bigl|\gamma(\beta_{xj})( {\mathrm e}^{-w_j^+} - 1 )\bigr|^2} \quad\text{and}\quad |c_j^-|^2 = \frac{4 |B|^2}{s_1^2\bigl|\gamma(\beta_{xj})( {\mathrm e}^{w_j^-} - 1 )\bigr|^2}, \quad j \in\mathcal{M}.$$ (51) On the other hand, if $$\gamma(\beta_{xj})$$ is real, then $$(w_j^-)^* = -w_j^+$$, so   $$\label{cjpcjm} {\text{Im}}\bigl[ c_j^+ (c_j^-)^*\bigr] = \left|\frac{B}{s_1 \gamma(\beta_{xj})}\right|^2 {\text{Im}} \bigl[{\text{csch}}^2( w_j^+/2)\bigr], \quad j \in\mathcal{N}.$$ (52) There is no dependence on $$q$$ here, which should be expected because Bloch waves propagate through the lattice without loss of energy. For a given value for $$\beta_y$$, we may now compute $$\langle E_{\mathcal S}\rangle$$ using (25), up to a factor of $$|B|^2$$, which plays no role in determining the direction of propagation. If $$\langle E_{\mathcal S} \rangle > 0$$, or $$\langle E_{\mathcal S} \rangle < 0$$, then the Bloch wave transports energy in the direction of increasing, or decreasing $$y$$, respectively. A Bloch wave with $$\langle E_{\mathcal S} \rangle=0$$ does not transport any energy in the $$y$$ direction. 7. The excitation problem We now consider a plane wave   $$u^{\mathrm i}({\mathbf{{r}}}) = {\mathrm e}^{{\mathrm i} k(x\cos\psi_0 +y\sin\psi_0)}$$ (53) incident on a semi-infinite lattice formed from small cylinders centred at the points $${\mathbf{{r}}}={\mathbf{{R}}}_{jp}$$, with $$j\in\mathbb{Z}$$ and $$p=0,1,\ldots$$ (see Fig. 1). The scattered field for this problem inherits the one-dimensional quasiperiodicity property (16) with $$\beta_x=k\cos\psi_0$$, and so we have   $$\label{def:us} u^{\mathrm s}({\mathbf{{r}}}) = \sum_{p=0}^\infty A_p \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} j s_1k\cos\psi_0} {H}_0^{(1)}(k|{\mathbf{{r}}}-{\mathbf{{R}}}_{jp}|),$$ (54) for as yet unknown coefficients $$A_p$$. It should be noted that the quasiperiodicity condition generally precludes the existence of surface waves along the edge at $$y=0$$, because $$0\in\mathcal{M}$$ for all angles of incidence (see (20)), meaning there is always at least one grating mode that transports energy into free space. Fig. 1. View largeDownload slide A plane wave incident on a semi-infinite lattice Fig. 1. View largeDownload slide A plane wave incident on a semi-infinite lattice We will apply the Dirichlet boundary condition (2) on the cylinders centred at the points $${\mathbf{{r}}}={\mathbf{{R}}}_{0q}$$, for $$q=0,1,\ldots$$ Boundary conditions elsewhere then follow by quasiperiodicity. The field incoming towards the cylinder located at $${\mathbf{{r}}}={\mathbf{{R}}}_{0q}$$, which consists of the incident wave and the radiation from all the other cylinders, is given by   $$\label{uiq} u^{\mathrm i}_q({\mathbf{{r}}}) = {\mathrm e}^{{\mathrm i} k(x\cos\psi_0 +y\sin\psi_0)} + \sum_{p=0}^\infty A_p \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} j s_1 k\cos\psi_0} {H}_0^{(1)}(k|{\mathbf{{r}}}-{\mathbf{{R}}}_{jp}|),$$ (55) where the prime indicates that the term in which $$j=0$$ and $$p=q$$ is omitted from the summation. The scattered response is   $$u^{\mathrm s}_q({\mathbf{{r}}}) = A_q{H}_0^{(1)}(k|{\mathbf{{r}}}-{\mathbf{{R}}}_{0q}|).$$ (56) Comparing these to (6) and (7), we find that $$A_q=-Z_0 u_q^{\mathrm i}({\mathbf{{r}}}_{0q})$$. Evaluating (55) at $${\mathbf{{r}}}={\mathbf{{r}}}_{0q}$$ then yields   $$\label{lin_sys1} A_q + Z_0 \sum_{p=0}^\infty A_p \sideset{}{'}\sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} j s_1 k\cos\psi_0} {H}_0^{(1)}(kR_{j,p-q}) = T_q, \quad q = 0,1,\ldots$$ (57) where   $$T_q = -Z_0 {\mathrm e}^{{\mathrm i} q k(\eta_1 \cos\psi_0 + \eta_2\sin\psi_0)}.$$ (58) If we now write   $$\label{def:Sp} S_p = \left\{\begin{array}{ll} \sigma_0(k\cos\psi_0) &\text{if } p = 0,\\[3pt] G_0(p{\mathbf{{s}}}_2;k\cos\psi_0) &\text{otherwise,} \end{array}\right.$$ (59) then (57) becomes   $$\label{lin_sys2} A_q + Z_0 \sum_{p=0}^\infty A_p S_{q-p} = T_q, \quad q = 0,1,\ldots$$ (60) In cases where no Bloch waves are excited, $$A_q\to0$$ as $$q\to\infty$$, and this system can be solved by truncation. No further conditions are needed. However, if Bloch waves are excited then $$\hat{A}_q\not\to 0$$, and (60) cannot be solved by truncation. Moreover, the solution is not unique, because we have yet to apply a radiation condition. The appropriate condition is to forbid Bloch modes that transport energy towards $$y=0$$ from appearing in the far field inside the lattice. This will be incorporated into the analysis as part of the filtering method, discussed in the next section. 8. The filtering method In a case where a single Bloch wave is excited, we can write   $$A_q = B {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}} + \hat{A}_q.$$ (61) The correct value for $$\beta_y$$ must yield a solution to (43) with $$\langle E_{\mathcal S}\rangle \ge 0$$ as described in section 6. This has two effects: it causes $$\hat{A}_q$$ to vanish as $$q\to\infty$$, and it ensures that the far field inside the lattice does not include any Bloch waves incoming from infinity, so that the radiation condition is satisfied. Now the decomposition (61) introduces a new unknown into the problem (the Bloch wave amplitude $$B$$), but this can be ‘filtered out’ by writing   $$A_q^{(1)} = \begin{cases} A_0 &\text{if } q = 0, \\ A_q - {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}}A_{q-1} &\text{otherwise}. \end{cases}$$ (62) For $$q>0$$, we then have $$A_q^{(1)} = \hat{A}_q - {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}}\hat{A}_{q-1}$$, which vanishes as $$q\to\infty$$. Solving the recurrence relation for $$A_q$$, we find that   $$A_q = \sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}}.$$ (63) There are now two ways in which we may obtain a system of equations for $$A_q^{(1)}$$. One possibility is to take phase-shifted differences between consecutive equations in (60). This was the preferred approach in (15), but the equation in which $$q=0$$ requires special treatment (because there is no preceding equation). If two Bloch waves are excited, the equations with $$q=0$$ and $$q=1$$ must be treated as special cases, and so on, meaning that the process is difficult to generalise to account for arbitrary numbers of Bloch waves. Therefore we use the alternative method (also suggested in (15)), which is to substitute (63) into (60). In this way, we obtain   $$\sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} + Z_0 \sum_{p=0}^\infty \sum_{j=0}^p A_j^{(1)} {\mathrm e}^{{\mathrm i} (p-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} S_{q-p} = T_q, \quad q=0,1,\ldots$$ (64) Again, this cannot be solved by truncation, because the inner sum in the second term involves $$A_0^{(1)},A_1^{(1)},\ldots$$ for every $$p$$. However, if we interchange the summations, we find that   $$\sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} + Z_0 \sum_{j=0}^\infty A_j^{(1)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}} \Gamma_{j-q}({\boldsymbol{\beta}}) = T_q,$$ (65) where   $$\Gamma_j({\boldsymbol{\beta}}) =\left\{ \begin{array}{ll} G_0^{(j,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}) &\text{if } j > 0,\\[3pt] G_0^{(1,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + \sigma_0(k\cos\psi_0) &\text{if } j = 0, \\[3pt] G_0^{(j,-1)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + G_0^{(1,\infty)}({\mathbf{{0}}};{\boldsymbol{\beta}}) + \sigma_0(k\cos\psi_0) &\text{if } j < 0. \end{array}\right.$$ (66) Since $$A_q^{(1)}\to0$$ as $$q\to\infty$$, equation (65) can be truncated at $$j=P$$ and $$q=P$$, say. The resulting finite system is then solved numerically, and (63) is used to compute the original coefficients. If $$P$$ is such that the discarded coefficients are numerically negligible, then the Bloch wave amplitude can be calculated using (61); thus   $$B = {\mathrm e}^{-{\mathrm i} P {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}} A_P.$$ (67) We now look to generalise the above procedure, to account for cases in which multiple Bloch waves are excited. In place of (61), we write   $$A_q = \hat{A}_q + \sum_{m=1}^\lambda B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(m)}}.$$ (68) Again, the correct choices for the $$y$$ components of the Bloch vectors are those values in the range (45) that yield solutions to (43) with $$\langle E_{\mathcal S}\rangle \ge 0$$, so that $$\hat{A}_q\to0$$ as $$q\to\infty$$, and the far field cannot include incoming Bloch waves. We also define   $$A_q^{(\lambda)} =\left\{ \begin{array}{ll} A_q & \text{if } q=0\text{ or } \lambda = 0, \\[3pt] A_q^{(\lambda-1)} - {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(\lambda)}} A_{q-1}^{(\lambda-1)} & \text{otherwise,} \end{array}\right.$$ (69) so that $$A_q^{(0)}$$ is the original unfiltered coefficient, whereas $$A_q^{(\lambda)}$$ has the first $$\lambda$$ Bloch waves filtered out for $$q\ge\lambda$$. Solving the recurrence relation for $$A_q^{(\lambda-1)}$$ yields   $$A_q^{(\lambda-1)} = \sum_{j=0}^q A_j^{(\lambda)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda)}},$$ (70) but we require an expression for $$A_q^{(0)}$$ in terms of $$A_q^{(\lambda)}$$, for substitution into the linear system (60). If we start with $$A_q^{(0)}$$ and use (70) twice, we obtain   \begin{align*} A_q^{(0)} &= \sum_{j=0}^q A_j^{(1)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} \\ &= \sum_{j=0}^q \sum_{p=0}^j A_p^{(2)} {\mathrm e}^{{\mathrm i}(j-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} \\ &= \sum_{p=0}^q A_p^{(2)} {\mathrm e}^{-{\mathrm i} p{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} \sum_{j=p}^q {\mathrm e}^{{\mathrm i} j {\mathbf{{s}}}_2\cdot({\boldsymbol{\beta}}^{(2)}-{\boldsymbol{\beta}}^{(1)})} \\ & = d_{2,1}\sum_{p=0}^q A_p^{(2)} {\mathrm e}^{{\mathrm i} (q-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} + d_{1,2}\sum_{p=0}^q A_p^{(2)} {\mathrm e}^{{\mathrm i}(q-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}}, \end{align*} where   $$d_{n,m} = \frac{1}{1-{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot({\boldsymbol{\beta}}^{(n)}-{\boldsymbol{\beta}}^{(m)})}}.$$ (71) Repeating this process, we find that   \begin{align*} A_n^{(0)} &= d_{2,1}\sum_{p=0}^q \sum_{j=0}^p A_j^{(3)} {\mathrm e}^{{\mathrm i}(p-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(3)}} {\mathrm e}^{{\mathrm i} (q-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} + d_{1,2}\sum_{p=0}^q \sum_{j=0}^p A_j^{(3)} {\mathrm e}^{{\mathrm i}(p-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(3)}} {\mathrm e}^{{\mathrm i}(q- p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} \\ &= d_{2,1}d_{3,1} \sum_{j=0}^q A_j^{(3)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(1)}} + d_{1,2} d_{3,2}\sum_{j=0}^q A_j^{(3)} {\mathrm e}^{{\mathrm i}(q-j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(2)}} \\ &\quad + d_{1,3}d_{2,3} \sum_{j=0}^q A_j^{(3)} {\mathrm e}^{{\mathrm i}(q- j){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(3)}}, \end{align*} since   \begin{equation*} d_{2,1}d_{1,3}+d_{1,2}d_{2,3} = d_{1,3}d_{2,3}. \end{equation*} Therefore we may postulate the general result   $$A_q^{(0)} = \sum_{m=1}^\lambda Q_m^\lambda \sum_{j=0}^q A_j^{(\lambda)} {\mathrm e}^{{\mathrm i} (q-j) {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}},$$ (72) where   $$Q_1^1 = 1 \quad\text{and}\quad Q_m^\lambda = \prod_{\substack{n=1\\ n\neq m}}^\lambda d_{n,m}, \quad \lambda > 1.$$ (73) This is correct in the cases where $$\lambda=1$$, $$2$$ and $$3$$. Let us attempt an inductive step by using (70) in (72). We find that   \begin{align*} A_q^{(0)} &= \sum_{m=1}^\lambda Q_m^\lambda \sum_{j=0}^q \sum_{p=0}^j A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(j-p){\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}} {\mathrm e}^{{\mathrm i} (q-j) {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}} \\ &= \sum_{m=1}^\lambda Q_m^\lambda {\mathrm e}^{{\mathrm i} q {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} \sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{-{\mathrm i} p {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}} \sum_{j=p}^q {\mathrm e}^{{\mathrm i} j{\mathbf{{s}}}_2\cdot({\boldsymbol{\beta}}^{(\lambda+1)}-{\boldsymbol{\beta}}^{(m)})} \\ &= \sum_{m=1}^\lambda Q_m^\lambda d_{\lambda+1,m}\sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} + \sum_{m=1}^\lambda Q_m^\lambda d_{m,\lambda+1} \sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}} \\ &= \sum_{m=1}^\lambda Q_m^{\lambda+1}\sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} + \biggl( \sum_{m=1}^\lambda Q_m^\lambda d_{m,\lambda+1} \biggr)\sum_{p=0}^q A_p^{(\lambda+1)} {\mathrm e}^{{\mathrm i}(q-p) {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(\lambda+1)}}. \end{align*} If we can show that   $$\sum_{m=1}^\lambda Q_m^\lambda d_{m,\lambda+1} = Q_{\lambda+1}^{\lambda+1},$$ (74) then (72) is established for all $$\lambda\in\mathbb{N}$$. Evidently (74) is true if $$\lambda=1$$. Otherwise it is equivalent to writing   $$\sum_{m=1}^\lambda \prod_{\substack{n=1 \\ n\neq m}}^\lambda \frac{ d_{n,m} }{ d_{n,\lambda+1}} = 1,$$ (75) and if we set $$z_m={\mathrm e}^{-{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}$$, the proof follows immediately on applying the result in the appendix. We may now substitute (72) into (60), recalling that $$A_q^{(0)}=A_q$$. In this way we find that the linear system for $$A_q^{(\lambda)}$$ is   $$\sum_{m=1}^\lambda Q_m^\lambda\left[ \sum_{j=0}^q A_j^{(\lambda)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(m)}} + Z_0 \sum_{j=0}^\infty A_j^{(\lambda)} {\mathrm e}^{{\mathrm i} (q-j){\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(m)}} \Gamma_{j-q}({\boldsymbol{\beta}}^{(m)})\right] = T_q, \quad q=0,1,\ldots$$ (76) This can be solved by truncation at $$j=P$$ and $$q=P$$, say, and the original coefficients can be computed using (72). We may also determine the Bloch wave amplitudes using (72). Thus, if $$A_j^{(\lambda)}$$ is negligible for $$j\ge P$$, then   $$A_q^{(0)} = \sum_{m=1}^\lambda \biggl( Q_m^\lambda \sum_{j=0}^P A_j^{(\lambda)} {\mathrm e}^{-{\mathrm i} j {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}}\biggr) {\mathrm e}^{{\mathrm i} q {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}}, \quad q \ge P,$$ (77) from which we can read off the Bloch wave amplitudes as   $$B^{(m)} = Q_m^\lambda \sum_{j=0}^P A_j^{(\lambda)} {\mathrm e}^{-{\mathrm i} j {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(m)}}}.$$ (78) In the case where $$\lambda=1$$ (so that $$m=1$$ is the only possibility), using (73) and (70) reduces this to   $$B^{(1)} = \sum_{j=0}^P A_j^{(\lambda)} {\mathrm e}^{-{\mathrm i} j {\mathbf{{s}}}_2 \cdot{{\boldsymbol{\beta}}^{(1)}}} = {\mathrm e}^{-{\mathrm i} P {\mathbf{{s}}}_2 \cdot {\boldsymbol{\beta}}^{(1)}} A_P^{(0)},$$ (79) in agreement with (67). 9. Infinite array subtraction Another idea from (15) offers a useful means of testing our results. If the Bloch vectors $${\boldsymbol{\beta}}^{(m)}$$ and associated amplitude coefficients $$B^{(m)}$$ have been determined, then we may form a linear system of equations for $$\hat{A}_q$$ by substituting (68) into (60), and taking known terms to the right-hand side. In this way, we obtain   $$\hat{A}_q + Z_0 \sum_{p=0}^\infty \hat{A}_p S_{q-p} = T_q - \sum_{m=1}^\lambda B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)} } \Bigl[ 1 + Z_0 \Gamma_{-q}\bigl({\boldsymbol{\beta}}^{(m)}\bigr)\Bigr], \quad q = 0,1,\ldots$$ (80) where $$\Gamma$$ is defined in (66). In this way, the slowly decaying terms are subtracted from the linear system, leaving a rapidly convergent system for the coefficients $$\hat{A}_q$$. The new system can be simplified slightly by exploiting the fact that each Bloch wave is a solution to the propagation problem. Thus, the dispersion relation (40) can be written in the form   $$1 + Z_0 \bigl[ G_0^{(-\infty,-q-1)}(0;{\boldsymbol{\beta}}) + \Gamma_{-q}({\boldsymbol{\beta}}) \bigr] = 0,$$ (81) so that   $$\hat{A}_q + Z_0 \sum_{p=0}^\infty \hat{A}_p S_{q-p} = T_q + Z_0 \sum_{m=1}^\lambda B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} G_0^{(-\infty,-q-1)}\bigl({\mathbf{{0}}};{\boldsymbol{\beta}}^{(m)}\bigr), \quad q = 0,1,\ldots$$ (82) It should be emphasised that this system will only converge if the correct values for $${\boldsymbol{\beta}}^{(m)}$$ and $$B^{(m)}$$ are inserted on the right-hand side. 10. Reflection and transmission It is convenient to represent the reflected and transmitted fields in terms of scattering angles $$\psi_j$$, which are defined so that   $$k\cos\psi_j = k\cos\psi_0 + 2j\pi/s_1 \quad\text{and}\quad k\sin\psi_j = {\mathrm i}\gamma(k\cos\psi_j),$$ (83) where $$\gamma$$ is given by (19). We also introduce the quantities   $$\rho_j = {\mathrm e}^{-{\mathrm i} k( \eta_1\cos\psi_j+\eta_2\sin\psi_j )} \quad\text{and}\quad \tau_j = {\mathrm e}^{{\mathrm i} k(\eta_1\cos\psi_j-\eta_2 \sin\psi_j)}.$$ (84) Note that $$\cos\psi_j$$ is always real, whereas $$\sin\psi_j$$ may be positive real or positive imaginary (but not zero since we are not dealing with Wood’s anomalies; see section 4). Consequently, $$|\rho_j|=|\tau_j|\to\infty$$ as $$|j|\to\infty$$. Using this notation, we can simplify the exponentials involving $$w_j^\pm$$, which appear in the spectral representations for the quasiperiodic Green’s functions (34) and (35). Thus, since $$\beta_x=k\cos\psi_0$$ for all Bloch vectors $${\boldsymbol{\beta}}$$, we find that   $${\mathrm e}^{w_j^+} = \rho_j {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}} \quad\text{and}\quad {\mathrm e}^{w_j^-} = \tau_j^{-1} {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}}.$$ (85) Below the lattice, the total field can be expanded as a series of grating modes as in (17), with $$\beta_{xj}=k\cos\psi_j$$. The only upward-oriented mode in this region is the incident field, so   $$u({\mathbf{{r}}}) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} kx\cos\psi_j} \bigl[ c_{j0}^- {\mathrm e}^{-{\mathrm i} ky\sin\psi_j} + \delta_{j0} {\mathrm e}^{{\mathrm i} ky\sin\psi_j} \bigr], \quad y \le 0,$$ (86) where $$\delta_{j0}=1$$ if $$j=0$$ and $$0$$ otherwise. To determine $$c_{j0}^-$$, we insert the quasiperiodic Green’s function (26) into the representation of the scattered field (54), to obtain   $$u^{\mathrm s}({\mathbf{{r}}}) = \sum_{p=0}^\infty A_p G_0({\mathbf{{r}}}-p{\mathbf{{s}}}_2;k\cos\psi_0).$$ (87) Then, introducing the decomposition (68) and using (31), we find that   $$u^{\mathrm s}({\mathbf{{r}}}) = \sum_{m=1}^\lambda B^{(m)} G_0^{(0,\infty)}\bigl({\mathbf{{r}}};{\boldsymbol{\beta}}^{(m)}\bigr) + \sum_{p=0}^\infty \hat{A}_p G_0({\mathbf{{r}}}-p{\mathbf{{s}}}_2;\beta_x).$$ (88) Finally, after using the spectral representations (27) and (34) along with (83) and (85), we may read off the reflection coefficients as   $$c_{j0}^- = {}\frac{2}{ks_1 \sin\psi_j}\Biggl[ \sum_{m=1}^\lambda \frac{B^{(m)}\tau_j }{ \tau_j-{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}} + \sum_{p=0}^\infty \hat{A}_p \tau_j^{-p} \Biggr].$$ (89) Evidently, the far field below the lattice may be obtained by simply limiting the summation in (86) to $$j\in\mathcal M$$, so that $$\sin\psi_j$$ is real, meaning that propagating modes are retained, whereas evanescent modes are discarded. Calculation of the transmitted field is more complicated. Between each consecutive pair of rows, the total field can again be represented in the form (17), this time with a full set of upward- and downward-oriented modes. That is,   $$u({\mathbf{{r}}}) = \sum_{j=-\infty}^\infty {\mathrm e}^{{\mathrm i} kx\cos\psi_j} \bigl[ c_{jq}^- {\mathrm e}^{-{\mathrm i} ky\sin\psi_j} + c_{jq}^+ {\mathrm e}^{{\mathrm i} ky\sin\psi_j} \bigr], \quad (q-1)\eta_2 \le y \le q\eta_2, \quad q \in\mathbb{N}.$$ (90) Since the far field inside the lattice can consist only of Bloch waves, we decompose the total field into the form   $$u({\mathbf{{r}}}) = \hat{u}({\mathbf{{r}}}) + u^{\mathrm{b}}({\mathbf{{r}}}),$$ (91) with   \begin{align} \hat{u}({\mathbf{{r}}}) = {\mathrm e}^{{\mathrm i} k(x\cos\psi_0+y\sin\psi_0)} - \sum_{m=1}^{\lambda} B^{(m)} G_0^{(-\infty,-1)}\bigl({\mathbf{{r}}};{\boldsymbol{\beta}}^{(m)}\bigr) + \sum_{p=0}^\infty \hat{A}_p G_0({\mathbf{{r}}}-p{\mathbf{{s}}}_2;\beta_x) \end{align} (92) and   \begin{align} u^{\mathrm{b}}({\mathbf{{r}}}) = \sum_{m=1}^{\lambda} B^{(m)}G_0^{(-\infty,\infty)}\bigl({\mathbf{{r}}};{\boldsymbol{\beta}}^{(m)}\bigr). \end{align} (93) We then write   $$c_{jq}^\pm = \hat{c}_{jq}^\pm + b_{jq}^\pm,$$ (94) where the terms on the right-hand side are the contributions from (92) and (93), respectively. For $$\hat{c}_{jq}^-$$, we observe that only the last term on the right-hand side of (92) includes downward-oriented modes (and only for $$p\ge q$$). Therefore, on using (27), we find that   $$\hat{c}_{jq}^- = \frac{2\tau_j^{-q}}{ks_1\sin\psi_j}\sum_{p=0}^\infty \hat{A}_{p+q} \tau_j^{-p}.$$ (95) One the other hand, all three terms in (92) contribute to $$\hat{c}_j^+$$. Again using (27), and also (35), we obtain   $$\hat{c}_{jq}^+ = \delta_{j0} + \frac{2}{ks_1\sin\psi_j}\biggl[ \sum_{m=1}^\lambda \frac{B^{(m)}} {1-\rho_j {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}} + \sum_{p=0}^{q-1} \hat{A}_p \rho_j^p \biggr].$$ (96) Similarly, for the contributions to $$u^\mathrm{b}$$, we use (34) with $$q_0=q$$ and (35) with $$q_1=q-1$$ to show that   $$b_{jq}^{+} = \frac{2\rho_j^q}{ks_1\sin\psi_j} \sum_{m=1}^\lambda \frac{B^{(m)} {\mathrm e}^{{\mathrm i} q{\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}{\rho_j{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}-1} \quad\text{and}\quad b_{jq}^{-} = \frac{2\tau_j^{-q}}{ks_1\sin\psi_j}\sum_{m=1}^\lambda \frac{ B^{(m)}{\mathrm e}^{{\mathrm i} q {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}} } {1-\tau_j^{-1}{\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}.$$ (97) We now consider the various contributions to the far field inside the lattice. It should be noted that, because each instance of (90) (that is, with a particular value for $$q$$) is valid for a finite range of $$y$$ values, we cannot simply discard evanescent modes. The behaviour of the amplitude coefficients as $$q\to\infty$$ must also be taken into account. Let us begin with (95). We write   $$y = (q-1) \eta_2 + y_0, \quad 0 \le y_0 \le \eta_2,$$ (98) so that the modes have the form   $$\hat{c}_{jq}^- {\mathrm e}^{{\mathrm i} k (x\cos\psi_j-y\sin\psi_j)} = \frac{2{\mathrm e}^{-{\mathrm i} q k \eta_1\cos\psi_j} }{ks_1\sin\psi_j}\sum_{p=0}^\infty \hat{A}_{p+q} \tau_j^{-p} {\mathrm e}^{{\mathrm i} k (x\cos\psi_j+(\eta_2-y_0)\sin\psi_j)}.$$ (99) This disappears as $$q\to\infty$$, due to the decay in the coefficients $$\hat{A}_{p+q}$$. Next consider (96), in the case where $$j\in\mathcal N$$. Since $$0\in\mathcal{M}$$, the first term on the right-hand side disappears. The second term has no dependence on $$q$$, and makes no contribution to the far field because it is to be multiplied by $${\mathrm e}^{-{\mathrm i} k y \sin\psi_j}$$, which decays exponentially as $$y$$ is increased. For the third term, we can proceed as above using (98), and write   $$\sum_{p=0}^{q-1} \hat{A}_p \rho_j^p \, {\mathrm e}^{{\mathrm i} k(x\cos\psi_j+y\sin\psi_j)} = {\mathrm e}^{{\mathrm i} k(x-(q-1)\eta_1)\cos\psi_j} {\mathrm e}^{{\mathrm i} k y_0\sin\psi_j} \sum_{p=0}^{q-1} \hat{A}_p \rho_j^{p-q+1}.$$ (100) In the far field, these contributions may be small because they originate from some distance away and decay exponentially in $$y$$, or, if they originate from nearby rows, because the amplitude coefficient $$\hat{A}_p$$ multiplying the quasiperiodic Green’s function is negligible. To see this mathematically, we split the series after roughly half the required terms. Thus,   $$\sum_{p=0}^{q-1} \hat{A}_p \rho_j^{p-q+1} = \sum_{p=0}^{{[(q-1)/2]}} \hat{A}_p \rho_j^{p-q+1} + \sum_{{p=[(q-1)/2]+1}}^{q-1} \hat{A}_p \rho_j^{p-q+1},$$ (101) where $$[\cdot]$$ means integer part (round towards $$-\infty$$). The first half of the sum disappears as $$q\to\infty$$ because $$|\rho_j|>1$$, and the second half disappears because $$|\rho_j|^{p-q+1}\le 1$$ and $$\hat{A}_p\to0$$ as $$p\to\infty$$. Finally, consider (96) for propagating modes. Since the far field consists of Bloch waves only, it must be the case that   $$\delta_{j0} + \frac{2}{ks_1\sin\psi_j}\biggl[ \sum_{m=1}^\lambda \frac{B^{(m)}} {1-\rho_j {\mathrm e}^{{\mathrm i} {\mathbf{{s}}}_2\cdot{\boldsymbol{\beta}}^{(m)}}} + \sum_{p=0}^\infty \hat{A}_p \rho_j^p \biggr] = 0, \quad j\in\mathcal{M}.$$ (102) In (5), the fact that these terms vanish identically arises directly from the analysis. Despite considerable effort on the part of the authors, we have not been able to prove this via the current method. Instead, we must view (102) as a necessary but not sufficient condition for the correctness of our numerical results, much like the conservation of energy condition derived in the next section. Ultimately, the far field as $$y\to\infty$$ is given by (90), with $$c_{jq}^\pm$$ replaced by $$b_{jq}^\pm$$ (that is, discarding $$\hat{c}_{jq}^\pm$$ from (94)). 11. Conservation of energy A conservation of energy condition for this problem was derived in (5), and we need not repeat the details here. The procedure is to evaluate the integral (21) around the parallelogram with vertices at   $${\mathbf{r}} = {}\pm \tfrac{1}{2}{\mathbf{s}}_1 \pm (N+\tfrac{1}{2}) {\mathbf{s}}_2, \quad N\in\mathbb{N},$$ (103) taking the derivative in the direction of the outgoing normal. The contributions from the sides parallel to $${\mathbf{s}}_2$$ cancel each other, due to the one-dimensional quasiperiodicity property (16). For the edge beneath the lattice, we simply evaluate (25) for the reflected field (86), taking into account the fact that the outgoing normal is $$-\hat{{\mathbf{y}}}$$, so that a factor of $$-1$$ is introduced. Since $$0\in\mathcal M$$ for any angle of incidence, the sum over $$\mathcal{N}$$ disappears, and we find that the time averaged energy flux across this line is   $$\label{E1} E_1 = {}-\frac{P_0\omega s_1 k}{2} \biggl[ \sin\psi_0 - \sum_{j\in\mathcal M} \sin\psi_j |c_{j0}^-|^2 \biggr],$$ (104) where $$c_{j0}^-$$ is given by (89). If no Bloch waves are excited, then all of the energy incident on the lattice must be reflected back, so we must have $$E_1=0$$. If Bloch waves are excited, then the average energy flux is in direction of the inward normal (some incident wave energy is transmitted into the lattice, so less is reflected back). In this case $$E_1<0$$ and   $$\label{cons_energy} E_1 + E_2 = 0,$$ (105) where $$E_2$$ is the time averaged energy flux across the upper edge of the parallelogram. To determine $$E_2$$, we again use (25), this time with the outgoing normal directed upward. We may assume that $$N$$ is sufficiently large to allow contributions from $$\hat{c}_{jq}^+$$ to be neglected (see section 10, above), so that the amplitude coefficients are given by (97). Some simplification then occurs because $$|\rho_j|=|\tau_j|=1$$ if $$j\in\mathcal{M}$$, whereas if $$j\in\mathcal{N}$$ then $$\tau_j^*=\rho_j$$. The average energy flux across the upper edge is therefore given by   \begin{align}\label{E2} E_2 &= \frac{2P_0\omega }{k s_1} \sum_{j\in\mathcal M} \frac{1}{\sin\psi_j} \left\{ \left| \sum_{m=1}^\lambda \frac{B^{(m)} {\mathrm e}^{{\mathrm i} N{\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}{\rho_j{\mathrm e}^{{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}-1}\right|^2 - \left|\sum_{m=1}^\lambda \frac{ B^{(m)}{\mathrm e}^{{\mathrm i} N{\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}} } {1-\tau_j^{-1}{\mathrm e}^{{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}\right|^2 \right\} \nonumber\\ &\quad - \frac{4P_0\omega}{k s_1} \sum_{j\in\mathcal N}\frac{1}{|\sin\psi_j|}\Im\left\{ \sum_{m=1}^\lambda \frac{B^{(m)} {\mathrm e}^{{\mathrm i} N{\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}{\rho_j{\mathrm e}^{{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}-1} \sum_{m=1}^\lambda \frac{ (B^{(m)})^*{\mathrm e}^{-{\mathrm i} N {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}} } {1-\rho_j^{-1}{\mathrm e}^{-{\mathrm i} {\mathbf{s}}_2\cdot{\boldsymbol{\beta}}^{(m)}}}\right\}. \end{align} (106) If a single Bloch wave is excited, then $$\lambda=1$$ and the factors involving $$N$$ all cancel. The meaning of this is important: we assumed $$N$$ to be large so that terms that do not contribute to the far field could be neglected. The resulting expression then consists only of contributions from Bloch waves. It does not depend on $$N$$ because Bloch waves propagate through the lattice without loss of energy. If multiple Bloch waves are excited then $$\lambda>1$$, and the explicit cancellation no longer occurs. However, $$E_2$$ still consists entirely of contributions from Bloch waves, and cannot depend on $$N$$. Therefore (106) can always be simplified by setting $$N=0$$, a fact which we have confirmed numerically. 12. Numerical results Figures 2 and 3 show the magnitudes of the coefficients $$A_j$$, $$A_j^{(\lambda)}$$ and $$\hat{A}_j$$ for parameters that lead to different behaviours as $$j$$ is increased. All of the results were computed using double precision arithmetic ($${\sim}16$$ decimal digits). For cases where Bloch waves are excited, a stringent test on accuracy is performed by computing $$A_j$$ using (76) and (72), and again by infinite array subtraction using (82) and (68), retaining only the Bloch wave amplitudes $$B^{(m)}$$ and vectors $${\boldsymbol{\beta}}^{(m)}$$ from the earlier calculation. We also confirm that (105) and (102) are satisfied in each case. Some care must be taken in assessing the accuracy of these, since the exact value is zero in both cases, so that the relative error is undefined. For (105), we use the number of cancelling digits in $$E_1$$ and $$E_2$$, or in cases where no Bloch waves are excited, the number of cancelling digits in the two terms in (104). For (102), we use the fact that there exists a propagating mode of unit amplitude (the incident field), so the number of correct significant figures is one greater than the number of correct decimal places. Fig. 2 View largeDownload slide Magnitudes of coefficients for a rectangular lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0,1]$$, $$a=0.005$$, $$\psi_0=0.25\pi$$ and (a) $$k=1.5$$, (b) $$k=3.0$$ Fig. 2 View largeDownload slide Magnitudes of coefficients for a rectangular lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0,1]$$, $$a=0.005$$, $$\psi_0=0.25\pi$$ and (a) $$k=1.5$$, (b) $$k=3.0$$ Fig. 3 View largeDownload slide Magnitudes of coefficients for a skewed lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0.1,1.2]$$, $$a=0.005$$, and (a) $$\psi_0=0.15\pi$$, $$k=3.7$$, (b) $$\psi_0=0.16\pi$$, $$k=3.525$$ Fig. 3 View largeDownload slide Magnitudes of coefficients for a skewed lattice with $${\mathbf{s}}_1=[1,0]$$, $${\mathbf{s}}_2=[0.1,1.2]$$, $$a=0.005$$, and (a) $$\psi_0=0.15\pi$$, $$k=3.7$$, (b) $$\psi_0=0.16\pi$$, $$k=3.525$$ Figure 2(a) shows $$|A_j|$$ in a case where no Bloch waves are excited. These results were obtained by solving (60) by truncation at $$p=200$$ and $$q=200$$, and may be used as a benchmark for the effectiveness of the filtering method. The rate of decay in the coefficients is such that the error introduced by truncation is extremely small; in fact $$|A_j|$$ is numerically negligible for $$j\gtrapprox45$$. The conservation of energy criterion (105), and (102) with $$j=0$$ (the only propagating mode), are both satisfied to near machine precision. In Fig. 2(b), the frequency has been increased so that a single Bloch wave is excited, with $$\beta_y\approx1.78$$. No other parameters have been changed. Evidently, solving (60) directly by truncation would fail, since $$A_j\not\to0$$ as $$j\to\infty$$. However, the filtered coefficients $$A_j^{(1)}\to 0$$ extremely rapidly, and are numerically negligible for $$j\gtrapprox10$$. Computing $$A_j$$ again using infinite array subtraction produces results in agreement to $$13$$–$$14$$ significant figures. There is a slight drop in accuracy for the conservation of energy criterion, in that now only $$13$$ digits in $$E_1$$ cancel those in $$E_2$$. Similarly, (102) is now satisfied to $$14$$ significant figures. This is to be expected, in view of the more complex calculations that are needed when Bloch waves are excited. Figure 3 shows results for two more difficult cases. In Fig. 3(a), the lattice is skewed, and two Bloch waves are excited, with $$\beta_y\approx0.919$$, and $$\beta_y\approx1.67$$. Again, the filtered coefficients $$A_j^{(2)}$$ and the coefficients $$\hat{A}_j$$ decay very rapidly as $$j$$ is increased, but only to a magnitude around $$10^{-12}$$, so that the results are slightly less accurate in this case. This is confirmed by computing $$A_j$$ again, using infinite array subtraction; the results were found to be in agreement to $$12$$ significant figures. For these parameters, there are two propagating modes ($$j=0$$ and $$j=1$$), and (102) was satisfied to $$12$$ significant figures in both instances, with similar accuracy in the conservation of energy criterion (105). Finally, in Fig. 3(b), the same lattice is used, but the wavenumber and angle of incidence are chosen so that a single Bloch wave is excited with $$\beta_y\approx1.48$$, and a second is close to cutting on (there would be two Bloch waves in the equivalent plot for $$k=3.526$$). A Bloch wave close to cutting on corresponds to a value for $$\beta_y$$ with a small, positive imaginary part. In this case, terms up to $$j=500$$ and $$p=500$$ were retained in the calculations because the decay in $$A_j^{(1)}$$ and $$\hat{A}_j$$ is very slow. In fact this ceases after $${\sim}250$$ rows, when $$|A_j^{(1)}|\approx|\hat{A}_j|\approx10^{-12}$$. For some rows, the agreement between $$A_j$$ computed by filtering and the value obtained by infinite array subtraction is as low as $$10$$ significant figures, but (102) (with $$j=0$$ and $$j=1$$) and the conservation of energy criterion (105) are still satisfied to $$12$$ significant figures. 13. Concluding remarks We have demonstrated a new method for calculating the reflection and transmission properties of lattices, without the need for integral transforms or the Wiener–Hopf technique, and with no spurious reflection or damping effects. The method yields very accurate results, with slight deterioration close to band edges and when multiple Bloch waves are excited. The former issue could be countered by including Bloch vectors with small imaginary parts in the filtering process (see section 8). Several of the steps undertaken in setting up the new method are also needed if the same problem is to be solved via the Wiener–Hopf technique. In particular, the starting point is the reduction of a boundary value problem to an infinite linear system of equations using multipoles, and the same lattice sums and quasiperiodic Green’s functions arise automatically. Solving the propagation problem to determine the Bloch vectors for waves that can be excited inside the lattice is necessary for both approaches. For the Wiener–Hopf technique, these Bloch vectors correspond to zeroes of the kernel (or zeroes of the determinant of the kernel) that lie on the unit circle. Information about these is needed in order to construct the correct inversion contour (see (5) for details). The Wiener–Hopf method also requires factorisation of the kernel function into a product of functions analytic inside and outside the unit circle. For point scatterer problems in which a single amplitude coefficient is associated with each lattice element, this can be performed by numerically evaluating a Cauchy integral as in (7), or by approximating the field between consecutive rows as a finite expansion in grating modes and then locating kernel zeroes inside the unit circle as in (5). On the other hand, our method requires only the numerical solution of a rapidly convergent infinite linear system. The Wiener–Hopf method has one advantage, in that the analysis proves that the far field inside the lattice is composed of Bloch waves only. We have been unable to prove that this is the case using the current method, but we have verified it numerically. In fact the requirement that all other contributions must vanish provides a useful check on our results. The differences between the two approaches are much more stark when finite size effects are considered. In this case, each lattice element is associated with a sequence of amplitude coefficients, and the quasiperiodic Green’s functions and lattice sums required are those considered in (6, Appendices A and B). The Wiener–Hopf equation for such problems has a matrix kernel which cannot be factorised by any known method. The approximate scheme used from (6) can be used instead, but this is complicated and very difficult to implement. On the other hand, our method generalises immediately. The algebra becomes more involved, but the procedure for obtaining solutions is largely identical to the point scatterer case. Semi-infinite lattices composed of finite-sized scatterers will be considered in the second part of this work. The technique can also be generalised to account for aperiodic sources (such as line sources) and lattices with defects, by using the array scanning method (22, 23), an approach which was developed for linear arrays in (24) and (25). References 1. Nicorovici N. A. McPhedran R. C. and Botten L. C. Photonic band gaps for arrays of perfectly conducting cylinders, Phys. Rev. E  52 ( 1995) 1135– 1145. Google Scholar CrossRef Search ADS   2. Sigalas M. and Economou E. N. Elastic and acoustic-wave band-structure, J. Sound Vib.  158 ( 1992) 377– 382. Google Scholar CrossRef Search ADS   3. Sigalas M. and Economou E. N. Band structure of elastic waves in two-dimensional systems, Solid State Commun.  86 ( 1993) 141– 143. Google Scholar CrossRef Search ADS   4. Joannopoulos J. D. Johnson S. G. Winn J. N. and Meade R. D. Photonic Crystals. Molding the Flow of Light , 2nd edn. ( Princeton University Press, Princeton, USA and Oxford, UK 2008). 5. Tymis N. and Thompson I. Low-frequency scattering by a semi-infinite lattice of cylinders, Q. J. Mech. Appl. Math.  64 ( 2011) 171– 195. Google Scholar CrossRef Search ADS   6. Tymis N. and Thompson I. Scattering by a semi-infinite lattice of cylinders and the excitation of Bloch waves, Q. J. Mech. Appl. Math.  67 ( 2014) 469– 503. Google Scholar CrossRef Search ADS   7. Haslinger S. G. Craster R. V. Movchan A. B. Movchan N. V. and Jones. I. S. Dynamic interfacial trapping of flexural waves in structured plates. Proc. R. Soc. A  472 ( 2016) 20150658 ( 25 pages). Google Scholar CrossRef Search ADS   8. Martin P. A. Multiple Scattering. Interaction of Time-Harmonic Waves with N Obstacles . ( Cambridge University Press, Cambridge, UK 2006). Google Scholar CrossRef Search ADS   9. Hills N. L. and Karp S. N. Semi-infinite diffraction gratings—I, Comm. Pure Appl. Maths  18 ( 1965) 203– 233. Google Scholar CrossRef Search ADS   10. Hills N. L. Semi-infinite diffraction gratings. II. Inward resonance, Comm. Pure Appl. Maths  18 ( 1965) 389– 395. Google Scholar CrossRef Search ADS   11. Linton C. M. and Martin P. A. Semi-infinite arrays of isotropic point scatterers. A unified approach, SIAM J. Appl. Math.  64 ( 2004) 1035– 1056. Google Scholar CrossRef Search ADS   12. Noble B. Methods Based on the Wiener–Hopf Technique . ( Chelsea, New York 1988). 13. Lawrie J. B. and Abrahams I. D. A brief historical perspective of the Wiener–Hopf technique. J. Engng. Math.  59 ( 2007) 351– 358. Google Scholar CrossRef Search ADS   14. Rogosin S. and Mishuris G. Constructive methods for factorization of matrix-functions, IMA J. Appl. Math.  81 ( 2016) 365– 391. Google Scholar CrossRef Search ADS   15. Linton C. M. Porter R. and Thompson I. Scattering by a semi-infinite periodic array and the excitation of surface waves, SIAM J. Appl. Math.  67, ( 2007) 1233– 1258. Google Scholar CrossRef Search ADS   16. Foldy L. L. The multiple scattering of waves I. General theory of isotropic scattering by randomly distributed scatterers, Phys. Rev.  67 ( 1945) 107– 119. Google Scholar CrossRef Search ADS   17. Brillouin L. Wave Propagation in Periodic Structures , 2nd edn. ( Dover, New York 1953). 18. Linton C. M. and Thompson I. Resonant effects in scattering by periodic arrays, Wave Motion  44 ( 2007) 167– 175. Google Scholar CrossRef Search ADS   19. Botten L. C. Nicorovici N. A. Asatryan A. A. McPhedran R. C. Martijn de Sterke C. and Robinson P. A. Formulation for electromagnetic scattering and propagation through grating stacks of metallic and dielectric cylinders for photonic crystal calculations. Part II. Properties and implementation. J. Opt. Soc. Am. A  17 ( 2000) 2177– 2190. Google Scholar CrossRef Search ADS   20. Twersky V. On the scattering of waves by an infinite grating, IRE Trans. on Antennas and Propagation  4 ( 1956) 330– 345. Google Scholar CrossRef Search ADS   21. Linton C. M. Lattice sums for the Helmholtz equation, SIAM Review  52 ( 2010) 630– 674. Google Scholar CrossRef Search ADS   22. Wu C. P. and Galindo V. Properties of a phased array of rectangular waveguides with thin walls, IEEE Trans. Antennas Propagat.  14 ( 1966) 163– 173. Google Scholar CrossRef Search ADS   23. Munk B. A. and Burrell G. A. Plane-wave expansion for arrays of arbitrarily oriented piecewise linear elements and its application in determining the impedance of a single linear antenna in a lossy half-space, IEEE Trans. Antennas Propagat.  27 ( 1979) 331– 343. Google Scholar CrossRef Search ADS   24. Thompson I. and Linton C. M. On the excitation of a closely spaced array by a line source, IMA J. Appl. Math.  72 ( 2007) 476– 497. Google Scholar CrossRef Search ADS   25. Thompson I. and Linton C. M. An interaction theory for scattering by defects in arrays, SIAM J. Appl. Math.  68 ( 2008) 1783– 1806. Google Scholar CrossRef Search ADS   Appendix Consider the function $$f$$, defined via   $$\label{def:f} f(z_1,\ldots,z_{\lambda+1}) = \sum_{m=1}^\lambda \prod_{\substack{n=1 \\ n \neq m}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_m}.$$ (A.1) We will prove that this is in fact the constant function $$f(z_1,\ldots,z_{\lambda+1})\equiv 1$$. First consider $$f$$ as a function of the complex variable $$z_1$$, with $$z_2,\ldots,z_{\lambda+1}$$ fixed (and distinct). For any integer $$p$$ between $$2$$ and $$\lambda$$ (inclusive), there are two ways in which $$z_1-z_p$$ can appear in the denominator: either $$m=1$$ or $$m=p$$. Separating these terms from the sum, and considering the limit in which $$z_1\to z_p$$, we find that   \begin{align*} \lim_{z_1\to z_p} \Biggl[ \prod_{n=2}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_1} + \prod_{\substack{n=1 \\ n \neq p}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_p} \Biggr] &= \prod_{\substack{n=2 \\ n \neq p}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_p} \lim_{z_1\to z_p} \Biggl[ \frac{z_p-z_{\lambda+1}}{z_p-z_1} + \frac{z_1-z_{\lambda+1}}{z_1-z_p} \Biggr] \\ &= \prod_{\substack{n=2 \\ n \neq p}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_p}. \end{align*} Therefore, $$f$$ is an entire function of the complex variable $$z_1$$. Furthermore, we may observe that   $$f(z_1,\ldots,z_{\lambda+1}) = \prod_{n=2}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_1} + \sum_{m=2}^\lambda \frac{z_1-z_{\lambda+1}}{z_1-z_m} \prod_{\substack{n=2 \\ n \neq m}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_m},$$ (A.2) from which we immediately find that   $$\lim_{z_1\to\infty } f(z_1,\ldots,z_{\lambda+1}) = \sum_{m=2}^\lambda \prod_{\substack{n=2 \\ n \neq m}}^\lambda \frac{z_n-z_{\lambda+1}}{z_n-z_m}.$$ (A.3) Hence, $$f$$ is a bounded entire function of $$z_1$$, so Liouville’s theorem shows that it must be independent of $$z_1$$. The same calculation can be made with $$z_n$$ in place of $$z_1$$, for $$n=2,\ldots,\lambda$$. Since replacing $$z_n$$ with $$z_n-\alpha$$ for $$n=1,\ldots,\lambda$$ does not affect $$f$$, we may make these translations in (A.1) to show that it is also independent of $$z_{\lambda+1}$$. Finally, we set $$z_1=z_{\lambda+1}$$, so that the terms with $$m>1$$ disappear from (A.1), and we easily find that   $$f(z_{\lambda+1},z_2,\ldots,z_{\lambda+1}) = 1,$$ (A.4) which completes the proof. © The Author, 2017. Published by Oxford University Press; all rights reserved. For Permissions, please email: journals.permissions@oup.com

### Journal

The Quarterly Journal of Mechanics and Applied MathematicsOxford University Press

Published: Feb 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations