# Diffraction of acoustic waves by a wedge of point scatterers111Submitted to
the SIAM Journal on Applied Mathematics on 06/08/2021.
Funding: A.V.K. is supported by a Royal Society Dorothy Hodgkin Research
Fellowship and a Dame Kathleen Ollerenshaw Fellowship, through which M.A.N.
has been funded.

M. A. Nethercote222Department of Mathematics, Alan Turing Building, University
of Manchester, United Kingdom,
<EMAIL_ADDRESS><EMAIL_ADDRESS>raphael.assier@manchester.ac.uk)., A. V. Kisil22footnotemark: 2 and R. C.
Assier22footnotemark: 2

###### Abstract

This article considers the problem of diffraction by a wedge consisting of two
semi-infinite periodic arrays of point scatterers. The solution is obtained in
terms of two coupled systems, each of which is solved using the discrete
Wiener–Hopf technique. An effective and accurate iterative numerical procedure
is developed to solve the diffraction problem, which allows us to compute the
interaction of thousands of scatterers forming the wedge. A crucial aspect of
this numerical procedure is the effective truncation of slowly convergent
single and double infinite series, which requires careful asymptotic analysis.
A convergence criteria is formulated and shown to be satisfied for a large
class of physically interesting cases. A comparison to direct numerical
simulations is made, highlighting the accuracy of the method.

##  1 Introduction

Wedge diffraction is a classical and well-studied problem Nethercote et al.
(2020). It has been considered in different media and with many different
boundary conditions (for hard, soft, impedance, penetrable and elastic
wedges). Diffraction by periodic structures is also an interesting and
challenging problem, which has been studied mainly for infinite Bruno and
Fernandez-Lado (2017); Movchan et al. (2018); Scarpetta and Sumbatyan (1996)
and semi-infinite arrays Hills and Karp (1965); Linton and Martin (2004). Much
less is known about diffraction by wedges made from periodic structures (in
particular point scatterers), which is the object of this article.

There has been renewed interest in periodic structures due to the growth of
the metamaterials research. Frequently in metamaterial research, some form of
periodicity is assumed which allows us to restrict the study to a small
portion known as the unit cell. The global behaviour is then reconstructed
from this local behaviour by symmetry. By modifying the unit cell, it is
possible to modify the properties of the whole material drastically. It can
display abnormal phenomena such as negative refraction, superlensing and
backward wave propagation which allow for the control and manipulation of
waves in ways that were originally perceived as unphysical. This has led to
many studies, which increase the complexity of the unit cell and reconstruct
the global scattering using the periodicity of the metamaterial Movchan et al.
(2018); Craster and Guenneau (2013).

An alternative approach looks into the case where periodicity is no longer
assumed. This means that the metamaterial will have well-defined boundaries
that can symbolise many different interfaces such as edges and corners. In
this framework, this article could be considered as diffraction by the
simplest possible wedge made out of a metasurface. Each face will be a semi-
infinite array made of sound-soft cylinders with infinite height and a small
radius, see Figure 2.

The wedge configuration can also be viewed as two separate semi-infinite
arrays with two sets of scattering coefficients to determine. From the two
derived infinite system of equations, we will find a solution using the
discrete Wiener–Hopf technique. The key steps of the technique are the
application of the Z-transform and the factorisation of the Wiener–Hopf kernel
(an infinite sum with some resemblance to the periodic Green’s function) on
the unit circle Camacho et al. (2019). Numerically, a number of challenges
need to be overcome due to slow converging series. Some tools that are
employed are tail-end asymptotics and rational approximations. It is found
that the latter approach captures the essential interesting behaviour and is
computationally fast. An iterative scheme is then constructed from the
individual Wiener–Hopf solutions.

One interesting feature of periodic structures is the possibility of resonance
effects. In the context of array scattering, it is frequently referred to as
Wood’s anomaly Liu and Declercq (2015). It was first observed in diffraction
gratings experimentally Wood (1902) and partly discussed theoretically by Lord
Rayleigh Rayleigh (1907). At resonance, the solution usually takes a simple
form of a plane wave moving along the array. Also, typically around the
resonant state, interesting phenomena can be observed such as, for example,
abnormal transmission or reflection Korolkov et al. (2016); Shanin and
Korolkov (2017). There are also different notions of resonance associated to
arrays (and no consistent terminology), some of which are the resonances of an
infinite array Evans_Porter_99 like Wood’s anomaly and some are associated
with effects of finiteness of the array Bennetts et al. (2017). A wedge with
periodic perturbations was considered in Kamotskiĭ and Nazarov (1999) and
Nazarov (2008). In this article, we will not be considering these cases since
they have distinctive features which require special attention and will be
considered in our future work.

The structure of the article is as follows. In Section 2, we start by
reviewing some fundamental canonical scattering problems related to infinite
and semi-infinite arrays. In particular, we review several important aspects
of multiple scattering theory, including Foldy’s approximation, isotropic
scattering and the discrete Wiener–Hopf technique. Section 3 focuses entirely
on diffraction by a wedge of point scatterers. After formulating the problem,
we derive an implicit solution via the discrete Wiener–Hopf technique and find
an approximation using iterative methods. This is followed by discussions on
the conditional convergence of the iterative scheme and comparisons with the
finite element software COMSOL MultiPhysics version 5.5. Lastly, in Section 4
we draw some conclusions and discuss further work. The appendix describes the
numerical methods, tail-end asymptotics and rational approximations used to
make the proposed methods fast and reliable computationally.

##  2 Scattering by infinite and semi-infinite arrays

In this article, we are investigating multiple scattering problems involving
semi-infinite arrays comprised of equidistant scatterers. The positions of the
scatterer centres are denoted by $\mathbf{R}_{n}(s,\alpha)$ where
$n\in{\mathbb{Z}}$ is the scatterer index, $s>0$ is the distance between
scatterers and $\alpha\in[0,\pi)$ is the angle between the array and the
$x$-axis. The explicit definition of $\mathbf{R}_{n}$ will change throughout
this article (see Figure 1 and 2). The scatterers themselves are all cylinders
of infinite height and small radius $a\in[0,s/2]$ and satisfy homogeneous
Dirichlet boundary conditions. The upper limit of $a$ is to prevent
overlapping cylinders. The incident wave field takes the form of a unit
amplitude plane wave at normal (non-skew) incidence. Problems of this type are
invariant with respect to the height position on the cylinder scatterers,
which simplifies the three-dimensional problem to a two-dimensional
projection.

We shall also be looking for time-harmonic solutions to the linear wave
equation by assuming and then suppressing the time factor $e^{-i\omega t}$
(where $\omega$ is the angular frequency). Let $\Phi$ be the total wave field,
which is decomposed into an incident field $\Phi_{\text{I}}$ and the resulting
scattered field $\Phi_{\text{S}}$ by the equation
$\Phi=\Phi_{\text{I}}+\Phi_{\text{S}}$ and use a polar coordinate system
$(r,\theta)$ with the position vector $\mathbf{r}$. The incident wave field
$\Phi_{\text{I}}(\mathbf{r})$ is given by,

$\displaystyle\Phi_{\text{I}}(\mathbf{r})=e^{-ikr\cos(\theta-\theta_{\text{I}})},$
(2.1)

where $k$ is the wavenumber and $\theta_{\text{I}}$ is the incoming incident
angle.

We assume that the cylinder scatterers are widely spaced out compared to their
size i.e. $a\ll s/2$ and we also assume that the cylinder radius is small in
comparison to the wavelength of the incident wave ($a\ll 2\pi/k$). These
assumptions alongside the Dirichlet boundary conditions allow us to make the
assumption of isotropic scattering in the spirit of Foldy (1945). This means
that for every cylinder, the scattered wave can be approximated by a Hankel
function of the first kind and zeroth order (denoted $H^{(1)}_{0}$) multiplied
by an unknown scattering coefficient. The total scattered field
$\Phi_{\text{S}}$ (which satisfies Helmholtz’s equation) is therefore a
summation of all these individual scattered waves,

$\displaystyle\Phi_{\text{S}}(\mathbf{r})=\sum_{n}\left[A_{n}H^{(1)}_{0}(k|\mathbf{r}-\mathbf{R}_{n}|)\right],$
(2.2)

where $A_{n}$ are the unknown scattering coefficients and the summation limits
will depend on the specifics of the array layout.

We can define the wave field incident on the $m^{\text{th}}$ scatterer,
$\Phi_{\text{I}}^{(m)}$, as the original incident wave plus the sum of all the
scattered waves from the rest of the cylinders. In other words, it is the
total wave field minus the $m^{\text{th}}$ scattered wave,

$\displaystyle\Phi_{\text{I}}^{(m)}(\mathbf{r})$
$\displaystyle=\Phi_{\text{I}}(\mathbf{r})+\sum_{n\neq
m}\left[A_{n}H^{(1)}_{0}(k|\mathbf{r}-\mathbf{R}_{n}|)\right],$
$\displaystyle=\Phi(\mathbf{r})-A_{m}H^{(1)}_{0}(k|\mathbf{r}-\mathbf{R}_{m}|).$
(2.3)

After applying the Dirichlet boundary conditions to (2), we use the assumption
that the cylinder radius $a$ is a small parameter to approximate,

$\displaystyle
A_{m}=-\frac{\Phi_{\text{I}}^{(m)}(\mathbf{R}_{m}+a\mathbf{\hat{r}})}{H^{(1)}_{0}(ka)}\approx-\frac{\Phi_{\text{I}}^{(m)}(\mathbf{R}_{m})}{H^{(1)}_{0}(ka)}.$
(2.4)

Note that normally $\Phi$ cannot be evaluated at $\mathbf{R}_{m}$ because it
diverges, but $\Phi_{\text{I}}^{(m)}$ can since the singularity has been
removed. We can use this to construct an infinite system of equations to solve
for the scattering coefficients,

$\displaystyle A_{m}H^{(1)}_{0}(ka)+\sum_{n\neq
m}\left[A_{n}H^{(1)}_{0}(k|\mathbf{R}_{m}-\mathbf{R}_{n}|)\right]+\Phi_{\text{I}}(\mathbf{R}_{m})=0.$
(2.5)

Before investigating the wedge problem in Section 3, we will review two well-
studied problems: the infinite array (see Figure 1(a)) and the semi-infinite
array (see Figure 1(b)). The scatterers of both problems have locations given
by,

$\displaystyle\mathbf{R}_{n}=ns(\cos(\alpha)\mathbf{\hat{x}}+\sin(\alpha)\mathbf{\hat{y}}),\
\text{where}\ \begin{cases}n=0,\pm 1,\pm 2,...&\text{infinite array,}\\\
n=0,1,2,...&\text{semi-infinite array,}\end{cases}$ (2.6)

in the Cartesian basis $(\mathbf{\hat{x}},\mathbf{\hat{y}})$.

Figure 1: Diagrams of the infinite (a) and semi-infinite (b) arrays with the
cylinders located at $\mathbf{R}_{n}$, the position vector $\mathbf{r}$ and
the incident wave $\Phi_{\text{I}}$.

For the infinite array problem (denoted with superscript ‘inf’), we have the
following infinite system of equations,

$\displaystyle
A^{\text{(inf)}}_{m}H^{(1)}_{0}\\!(ka)+\\!\\!\sum_{\begin{subarray}{c}n=-\infty\\\
n\neq
m\end{subarray}}^{\infty}\\!\left[A^{\text{(inf)}}_{n}H^{(1)}_{0}\\!(ks|m-n|)\right]\\!+e^{-iksm\cos(\theta_{\text{I}}-\alpha)}\\!=\\!0,\
\ m\in{\mathbb{Z}}.$ (2.7)

Since
$\Phi_{\text{I}}(\mathbf{r}+\mathbf{R}_{m})=e^{-iksm\cos(\theta_{\text{I}}-\alpha)}\Phi_{\text{I}}(\mathbf{r})$,
it is logical to say that $\Phi$ (and hence $A^{\text{(inf)}}_{m}$) satisfies
the same type of quasi-periodicity. We can therefore use
$A^{\text{(inf)}}_{m}=e^{-iksm\cos(\theta_{\text{I}}-\alpha)}A^{\text{(inf)}}_{0}$
so that (2.7) is reduced to solving for a single coefficient
$A^{\text{(inf)}}_{0}$ which is then given by,111The formula (2.8) is very
slow to converge, which is why, in practice, we sought an alternative formula
by using the same identities that lead to (A.7). The result is also quite
similar to equation (2.47) in Linton (1998)

$\displaystyle
A^{\text{(inf)}}_{0}=-\left(H^{(1)}_{0}(ka)+2\sum_{n=1}^{\infty}\left[\cos(ksn\cos(\theta_{\text{I}}-\alpha))H^{(1)}_{0}(ksn)\right]\right)^{-1}.$
(2.8)

Similarly, for the semi-infinite array problem (with superscript ‘s-inf’), we
have $A^{\text{(s-inf)}}_{m}$ satisfying the following infinite system of
equations,

$\displaystyle
A^{\text{(s-inf)}}_{m}H^{(1)}_{0}\\!(ka)+\\!\\!\sum_{\begin{subarray}{c}n=0\\\
n\neq
m\end{subarray}}^{\infty}\\!\left[A^{\text{(s-inf)}}_{n}H^{(1)}_{0}\\!(ks|m-n|)\right]\\!+e^{-iksm\cos(\theta_{\text{I}}-\alpha)}\\!=\\!0,\
\ m\geq 0,$ (2.9)

and $A^{\text{(s-inf)}}_{m}=0$ for all $m<0$.

We solve this system using the discrete analogue of the Wiener–Hopf technique.
Here, we will only give an overview of the technique for brevity (see Hills
and Karp (1965) and Linton and Martin (2004) for more details). We start by
extending the system of equations (2.9) for negative $m$ using some unknown
coefficients $F_{m}$,

$\displaystyle
A^{\text{(s-inf)}}_{m}H^{(1)}_{0}\\!(ka)+\\!\\!\sum_{\begin{subarray}{c}n=-\infty\\\
n\neq
m\end{subarray}}^{\infty}\\!\left[A^{\text{(s-inf)}}_{n}H^{(1)}_{0}\\!(ks|m-n|)\right]\\!=\\!\begin{cases}-e^{-iksm\cos(\theta_{\text{I}})}&m\geq
0,\\\ F_{m}&m<0.\end{cases}$ (2.10)

Then we multiply (2.10) by $z^{m}$ and summing from $m=-\infty$ to $\infty$
(i.e. the discrete $Z$-transform) to obtain the Wiener–Hopf equation,

$\displaystyle K(z)A^{+}(z)=F^{+}(z)+F^{-}(z),$ (2.11)

where the goal is to determine the $Z$-transform of the scattering
coefficients $A^{+}(z)$,

$\displaystyle
A^{+}(z)=\sum_{m=-\infty}^{\infty}A^{\text{(s-inf)}}_{m}z^{m}=\sum_{m=0}^{\infty}A^{\text{(s-inf)}}_{m}z^{m}.$
(2.12)

To help with the convergence of the $Z$-transform, it is useful to assume that
$k$ has a small positive imaginary part (i.e.
$0<\textup{Im}\left\\{k\right\\}\ll 1$). Note that we will let
$\textup{Im}\left\\{k\right\\}\rightarrow 0$ after we finish using the
Wiener–Hopf technique because the final solution is still valid when
$\textup{Im}\left\\{k\right\\}=0$. This will be consistent with all future
numerical computations presented in this article.

We use the superscripts ‘$+$’ (resp. ‘$-$’) for functions that are analytic
and zero-free in the region $\Omega^{+}$ (resp. $\Omega^{-}$), where,

$\displaystyle\Omega^{+}=\left\\{z\in{\mathbb{C}}:|z|<e^{-\textup{Im}\left\\{k\right\\}s\cos(\theta_{\textrm{I}}-\alpha)}\right\\},\
\
\Omega^{-}=\left\\{z\in{\mathbb{C}}:|z|>e^{-\textup{Im}\left\\{k\right\\}s}\right\\}.$
(2.13)

Note that the union $\Omega^{+}\cup\Omega^{-}$ is the entire complex $z$-plane
and the intersection $\Omega^{+}\cap\Omega^{-}$ is the annulus region
$e^{-\textup{Im}\left\\{k\right\\}s}<|z|<e^{-\textup{Im}\left\\{k\right\\}s\cos(\theta_{\textrm{I}}-\alpha)}$.
In the Wiener–Hopf equation (2.11), $F^{+}(z)$ is the discrete $Z$-transform
of the incident wave forcing and $F^{-}(z)$ is an unknown forcing term:

$\displaystyle
F^{+}(z)=\frac{1}{ze^{-iks\cos(\theta_{\text{I}}-\alpha)}-1},\quad
F^{-}(z)=\sum_{m=-\infty}^{-1}F_{m}z^{m},$ (2.14)

where the simple pole at $z=e^{iks\cos(\theta_{\text{I}}-\alpha)}$ is
considered to be outside the region $\Omega^{+}$. $K(z)$ is the discrete
Wiener–Hopf kernel given by,

$\displaystyle
K(z)=H^{(1)}_{0}(ka)+\sum_{\ell=1}^{\infty}\left[(z^{\ell}+z^{-\ell})H^{(1)}_{0}(ks\ell)\right].$
(2.15)

The Wiener–Hopf kernel has the property $K(z)=K(1/z)$ for all $z$ and is both
analytic and zero-free on the annulus region
$e^{-\textup{Im}\left\\{k\right\\}s}<|z|<e^{\textup{Im}\left\\{k\right\\}s}$
(which includes the intersection region $\Omega^{+}\cap\Omega^{-}$) except for
the branch point singularities at $z=e^{\pm iks}$. Note that (2.15) is not
well-suited for numerical evaluation. This is because the formula only
converges on
$e^{-\textup{Im}\left\\{k\right\\}s}<|z|<e^{\textup{Im}\left\\{k\right\\}s}$
and converges very slowly (see (A.2)). In Appendix A, we will discuss some
alternative methods to evaluate $K(z)$ including a different formula (A.7) and
an approximation (A.11) which will be essential for producing the results of
this article.

To proceed with the Wiener–Hopf technique, (2.11) needs to be rearranged such
that each side of the equality is both analytic and zero-free in $\Omega^{+}$
or $\Omega^{-}$. A key step of the technique is the factorisation of the
kernel, $K(z)=K^{+}(z)K^{-}(z)$, which is done such that the factors satisfy
the property $K^{+}(z)=K^{-}(1/z)$ for all $z$. This factorisation is achieved
by using Cauchy’s integral formula to sum split the natural logarithm of the
kernel

$\displaystyle\ln(K(z))$ $\displaystyle=\frac{1}{2\pi
i}\int_{C^{+}}\frac{\ln(K(\xi))}{\xi-z}\text{d}\xi-\frac{1}{2\pi
i}\int_{C^{-}}\frac{\ln(K(\xi))}{\xi-z}\text{d}\xi,$ (2.16)

where the integration contours $C^{\pm}$ are anticlockwise circular paths
contained inside $\Omega^{+}\cap\Omega^{-}$ on the $\xi$ complex plane.
Additionally $C^{+}$ (resp. $C^{-}$) will also run radially above (resp.
below) the pole at $\xi=z$. The former (resp. latter) integral is analytic in
$\Omega^{+}$ (resp. $\Omega^{-}$). If we apply a simple transformation
($\xi\rightarrowtail 1/\xi$) to the $C^{+}$ integral, we can rewrite (2.16)
such that the factors will obviously satisfy $K^{+}(z)=K^{-}(1/z)$,

$\displaystyle\ln(K^{+}(z))$ $\displaystyle=\ln(K_{0})-\frac{1}{2\pi
i}\int_{C^{-}}\frac{z\ln(K(\xi))}{z\xi-1}\text{d}\xi,$ (2.17)
$\displaystyle\ln(K^{-}(z))$ $\displaystyle=\ln(K_{0})-\frac{1}{2\pi
i}\int_{C^{-}}\frac{\ln(K(\xi))}{\xi-z}\text{d}\xi,$ (2.18)
$\displaystyle\text{where,}\quad\ln(K_{0})$
$\displaystyle=\ln(K^{+}(0))=\frac{1}{4\pi
i}\int_{C^{-}}\frac{\ln(K(\xi))}{\xi}\text{d}\xi.$ (2.19)

We use pole removal to sum split $\frac{F^{+}(z)}{K^{-}(z)}$ after factorising
the kernel and then apply Liouville’s theorem to solve the Wiener–Hopf
equation to obtain the transformed solution,

$\displaystyle
A^{+}(z)=\frac{1}{K^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})K^{+}(z)(ze^{-iks\cos(\theta_{\text{I}}-\alpha)}-1)}.$
(2.20)

Now, we apply the identity $K^{-}(z)=K^{+}(1/z)$ and the inverse
$Z$-transform,

$\displaystyle A^{\text{(s-inf)}}_{m}=\frac{1}{2\pi
i}\int_{C^{-}}A^{+}(z)z^{-m-1}\text{d}z,$ (2.21)

to determine the scattering coefficients for the semi-infinite array problem,

$\displaystyle
A^{\text{(s-inf)}}_{m}=-\frac{1}{K^{+}(e^{-iks\cos(\theta_{\text{I}}-\alpha)})}\sum_{n=0}^{m}\left[\lambda_{n}e^{-iks(m-n)\cos(\theta_{\text{I}}-\alpha)}\right].$
(2.22)

Here the coefficients $\lambda_{n}$ are defined by the following Taylor
series,

$\displaystyle\frac{1}{K^{+}(z)}=\sum_{n=0}^{\infty}\lambda_{n}z^{n},\ \
\text{where}\ \
\lambda_{n}=\frac{1}{n!}\frac{\text{d}^{n}}{\text{d}z^{n}}\left[\frac{1}{K^{+}(z)}\right]_{z=0},$
(2.23)

and can be given using Cauchy’s integral formula,

$\displaystyle\lambda_{n}=\frac{1}{2\pi
i}\int_{C^{-}}\frac{z^{-n-1}}{K^{+}(z)}\text{d}z=\frac{1}{2\pi}\int_{-ks}^{2\pi-
ks}\frac{e^{-int}}{K^{+}(e^{it})}\text{d}t.$ (2.24)

It should be noted that $\frac{1}{K^{+}(e^{it})}$ is $L^{1}$ integrable on the
real line segment $[-ks,2\pi-ks]$ (i.e. the integral
$\int_{0}^{2\pi}\left|\frac{1}{K^{+}(e^{it})}\right|\text{d}t$ is finite). By
the Riemann-Lebesgue lemma, we can therefore say that $\lambda_{n}\rightarrow
0$ as $n\rightarrow\infty$. The $\lambda_{n}$ coefficients for $n\geq 1$ can
also be written without requiring an exact formula for the $K^{+}$ factor. We
do this by using the original definition (2.23), the following identity,

$\displaystyle\frac{\text{d}}{\text{d}z}\left[\frac{1}{K^{+}(z)}\right]=-\frac{1}{K^{+}(z)}\frac{\text{d}}{\text{d}z}\left[\ln\left(K^{+}(z)\right)\right],$
(2.25)

and the derivatives of (2.17), which results in the following formula,

$\displaystyle\lambda_{n}=-\sum_{m=1}^{n}\frac{m\lambda_{n-m}}{n}\int_{0}^{1}\cos(m\pi\tau)\ln\left(K(e^{i\pi\tau})\right)\text{d}\tau,\quad
n\geq 1.$ (2.26)

Lastly, we can also rewrite the equation for the scattering coeffcients (2.22)
as a recurrence relation,

$\displaystyle A^{\text{(s-inf)}}_{m}$
$\displaystyle\\!=\\!e^{-iks\cos(\theta_{\text{I}}-\alpha)}A^{\text{(s-inf)}}_{m-1}\\!+\\!\frac{\lambda_{m}}{\lambda_{0}}A^{\text{(s-inf)}}_{0},\
\
A^{\text{(s-inf)}}_{0}\\!=\\!-\frac{\lambda_{0}}{K^{+}(e^{-iks\cos(\theta_{\text{I}}-\alpha)})}\cdot$
(2.27)

This recurrence relation is useful because it is more efficient for computing
the scattering coefficients $A^{\text{(s-inf)}}_{m}$ and it highlights their
asymptotic behaviour as $m\rightarrow\infty$. Specifically, they become
equivalent to the solution of the infinite array problem
$A^{\text{(s-inf)}}_{m}\sim A^{\text{(inf)}}_{m}$ as $m\rightarrow\infty$.

##  3 Diffraction by a wedge of point scatterers

Now that we have discussed the infinite and semi-infinite array problems, we
proceed to the wedge problem (see Figure 2). Here we consider each wedge face
to be formed from a semi-infinite array of cylinder scatterers. The positions
of these scatterers are given by,

$\displaystyle\mathbf{R}_{n}=$
$\displaystyle|n|\begin{cases}\mathbf{s}_{+},&n\geq 0,\\\
\mathbf{s}_{-},&n<0,\end{cases}=|n|s\begin{cases}\cos(\alpha)\mathbf{\hat{x}}+\sin(\alpha)\mathbf{\hat{y}},&n\geq
0,\\\
\cos(\alpha)\mathbf{\hat{x}}-\sin(\alpha)\mathbf{\hat{y}},&n<0,\end{cases}$
(3.1)

for $n\in{\mathbb{Z}}$. To prevent the cylinders overlapping, this
configuration requires $a<s/2$ (as previously stated in Section 2) as well as
the additional condition $\sin(\alpha)>a/s$. Aside from $\mathbf{R}_{n}$, we
will use the same definitions and assumptions as Section 2. For the total
field $\Phi$, we pose a solution of the form,

$\displaystyle\Phi(\mathbf{r})=\Phi_{\text{I}}+\sum_{n=0}^{\infty}\left[A_{n}H^{(1)}_{0}(k|\mathbf{r}-\mathbf{R}_{n}|)\right]+\sum_{n=-\infty}^{-1}\left[B_{n}H^{(1)}_{0}(k|\mathbf{r}-\mathbf{R}_{n}|)\right],$
(3.2)

where $A_{n}$ and $B_{n}$ are the coefficients for scatterers on the upper and
lower arrays $n\geq 0$ and $n<0$ respectively. Because of this restriction, we
set $A_{n}=0$ for all $n<0$ and $B_{n}=0$ when $n\geq 0$.

Figure 2: Diagram of the point scatterer wedge with scatterers located at
$\mathbf{R}_{n}$ (see (3.1)), the position vector $\mathbf{r}$ and the
incident wave $\Phi_{\text{I}}$.

Initially, we shall approximate the scattering coefficients $A_{n}$ and
$B_{n}$ by assuming that the two arrays are isolated from each other. This
gives us two separate semi-infinite array problems which we can relatively
easily find an exact solution for,

$\displaystyle A_{n}\approx A^{(0)}_{n}$
$\displaystyle=-\frac{e^{-iksn\cos(\theta_{\text{I}}-\alpha)}}{K^{+}(e^{-iks\cos(\theta_{\text{I}}-\alpha)})}\sum_{m=0}^{n}\lambda_{m}e^{iksm\cos(\theta_{\text{I}}-\alpha)}\equiv
A^{\text{(s-inf)}}_{n},$ (3.3) $\displaystyle B_{-n}\approx B^{(0)}_{-n}$
$\displaystyle=-\frac{e^{-iksn\cos(\theta_{\text{I}}+\alpha)}}{K^{+}(e^{-iks\cos(\theta_{\text{I}}+\alpha)})}\sum_{m=0}^{n-1}\lambda_{m}e^{iksm\cos(\theta_{\text{I}}+\alpha)}.$

The isolation approximation given by (3.3) allows us to get an early glimpse
into the behaviour of the scattering coefficients as $n$ becomes large and
will be quite useful later in this article as a preliminary approximation.
Figure 3 illustrates an example of this semi-infinite array approximation
including scattering coefficients up to and including $A_{1000}$ and
$B_{-1000}$.

Figure 3: The semi-infinite array approximation (3.3) with the wedge defined
by $\alpha=\frac{5\pi}{6}$, $s=0.1$ and $a=0.01$. Here we plot the real part
of the scattered wave with the left (resp. right) side having an incident wave
defined by $k=5\pi$ and $\theta_{\textrm{I}}=0$ (resp. $k=15\pi$ and
$\theta_{\textrm{I}}=\frac{\pi}{2}$).

Following the same reasoning as in Section 2, starting from (2.2) and applying
the boundary conditions, we can form two coupled infinite systems of equations
depending on the value of the index $m$ to solve for the scattering
coefficients $A_{m}$ and $B_{m}$. If the index is positive ($m\geq 0$), then
we solve a system of equations for the top array coefficient $A_{m}$ where the
forcing has dependency on the incident wave and the bottom array coefficients,

$\displaystyle A_{m}H^{(1)}_{0}(ka)$
$\displaystyle+\sum_{\begin{subarray}{c}n=0\\\ n\neq
m\end{subarray}}^{\infty}\left[A_{n}H^{(1)}_{0}(ks|m-n|)\right]$ (3.4)
$\displaystyle=-\sum_{n=-\infty}^{-1}\left[B_{n}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(m,-n)\right)\right]-e^{-iksm\cos(\theta_{\text{I}}-\alpha)},\quad
m\geq 0,$ $\displaystyle\textrm{where}\ \ \Lambda_{\alpha}(m$
$\displaystyle,n)=(m^{2}+n^{2}-2mn\cos(2\alpha))^{\frac{1}{2}}.$ (3.5)

If the index is strictly negative ($m<0$), then we solve for $B_{m}$ in terms
of the coefficients $A_{n}$ and $B_{n}$,

$\displaystyle B_{m}H^{(1)}_{0}(ka)$
$\displaystyle+\sum_{\begin{subarray}{c}n=-\infty\\\ n\neq
m\end{subarray}}^{-1}\left[B_{n}H^{(1)}_{0}(ks|m-n|)\right]$ (3.6)
$\displaystyle=-\sum_{n=0}^{\infty}\left[A_{n}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(-m,n)\right)\right]-e^{iksm\cos(\theta_{\text{I}}+\alpha)},\quad
m<0.$

From left to right, each of the terms in these two systems have their own
physical meaning: the scattered wave from the $m$th cylinder, the interaction
with the rest of the cylinders in the associated array, the forcing due to the
other array and the forcing due to the incident wave. Note that if we neglect
the third terms in these systems (the forcing due to the other array), we are
left with the semi-infinite array problems resulting from the isolation
assumption (see (3.3)). Solving these two systems will give us $A_{m}$ in
terms of $B_{m}$ and vice versa.

###  3.1 Discrete Wiener–Hopf technique on top array scattering coefficients

We can solve the top array system of equations (3.4) using the discrete
Wiener–Hopf technique in the same way as the semi-infinite array problem. To
start, we need to extend (3.4) using some unknown coefficients $F_{m}$ and
recalling $A_{m}=0$ for $m<0$,

$\displaystyle A_{m}H^{(1)}_{0}(ka)$
$\displaystyle+\sum_{\begin{subarray}{c}n=-\infty\\\ n\neq
m\end{subarray}}^{\infty}\\!\left[A_{n}H^{(1)}_{0}\\!(ks|m-n|)\right]$ (3.7)
$\displaystyle=\begin{dcases}-\sum_{n=-\infty}^{-1}\\!\left[B_{n}H^{(1)}_{0}\\!\left(ks\Lambda_{\alpha}(m,-n)\right)\right]\\!-e^{-iksm\cos(\theta_{\text{I}}-\alpha)},&m\geq
0,\\\ F_{m},&m<0.\end{dcases}$

Here the $A_{m}$ scattering coefficients (which are zero for all $m<0$, as
stated previously) are the unknowns to find and the $B_{m}$ coefficients
(which are zero for all $m\geq 0$) are assumed to be known. To obtain the
Wiener–Hopf equation, we apply the same Z-transform as in the semi-infinite
array problem,

$\displaystyle G(z)=\sum_{m=-\infty}^{\infty}G_{m}z^{m},\quad
G_{m}=\frac{1}{2\pi i}\oint_{|z|=1}G(z)z^{-m-1}\text{d}z,$ (3.8)

to (3.7) to obtain the Wiener–Hopf equation

$\displaystyle K(z)A^{+}(z)=$ $\displaystyle
F_{\text{pole}}^{+}+F_{B}^{+}(z)+F^{-}(z),$ (3.9)

where $A^{+}(z)$ is the transform of the unknown scattering coefficients given
by

$\displaystyle
A^{+}(z)=\sum_{m=-\infty}^{\infty}A_{m}z^{m}=\sum_{m=0}^{\infty}A_{m}z^{m}.$
(3.10)

Like the semi-infinite array problem, it is useful to assume that $k$ has a
small positive imaginary part to help with the convergence of the Z-transform.
We also reuse the definitions of $\Omega^{\pm}$ (see (2.13)), to indicate the
analytic regions of a function with a $+$ or $-$ superscript.

Of the other components in this Wiener–Hopf equation, the kernel $K(z)$ has
the exact same definition and properties as (2.15) including the important
identity $K(z)=K(1/z)$ and the singular points $z=e^{\pm iks}$. The three
forcing terms on the right-hand side are defined by

$\displaystyle F_{\text{pole}}^{+}(z)$
$\displaystyle=-\sum_{m=0}^{\infty}e^{-iksm\cos(\theta_{\text{I}}-\alpha)}z^{m}=\frac{1}{ze^{-iks\cos(\theta_{\text{I}}-\alpha)}-1},$
(3.11) $\displaystyle F_{B}^{+}(z)$
$\displaystyle=-\sum_{m=0}^{\infty}\sum_{n=1}^{\infty}B_{-n}z^{m}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(m,n)\right),$
(3.12) $\displaystyle F^{-}(z)$
$\displaystyle=\sum_{m=-\infty}^{-1}F_{m}z^{m}.$ (3.13)

Here, $F^{-}(z)=O\left(\frac{1}{z}\right)$ as $|z|\rightarrow\infty$ by
design.

We proceed with the Wiener–Hopf technique by factorising $K(z)$ in the exact
same way as Section 2. This means writing $K(z)=K^{+}(z)K^{-}(z)$ in such a
way that the two factors satisfy $K^{-}(1/z)=K^{+}(z)$ and are defined by
Cauchy’s integral formula (2.17)-(2.19). Now let’s divide (3.9) by $K^{-}(z)$,

$\displaystyle K^{+}(z)A^{+}(z)=$
$\displaystyle\frac{F_{\text{pole}}^{+}(z)}{K^{-}(z)}+\frac{F_{B}^{+}(z)}{K^{-}(z)}+\frac{F^{-}(z)}{K^{-}(z)}.$
(3.14)

Next, we sum-split each of the following functions,

$\frac{F_{\text{pole}}^{+}(z)}{K^{-}(z)}\ \ \text{and}\ \
\frac{F_{B}^{+}(z)}{K^{-}(z)},$

into two parts that are analytic on $\Omega^{\pm}$. For the first function,
sum-splitting is done by pole removal,

$\displaystyle\frac{F_{\text{pole}}^{+}(z)}{K^{-}(z)}=$ $\displaystyle\
\frac{F_{\text{pole}}^{+}(z)}{K^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})}+F_{\text{pole}}^{+}(z)\left[\frac{1}{K^{-}(z)}-\frac{1}{K^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})}\right].$
(3.15)

Here the first term of the right-hand side is analytic on $\Omega^{+}$ and the
second term is analytic on $\Omega^{-}$. The second function
$\frac{F_{B}^{+}(z)}{K^{-}(z)}$ is sum-split by rewriting it as a Laurent’s
series because it is analytic on a thin annulus region of $z$ containing
$\Omega^{+}\cap\Omega^{-}$.

$\displaystyle\frac{F_{B}^{+}(z)}{K^{-}(z)}=D(z)=\sum_{n=-\infty}^{\infty}D_{n}z^{n},\
\ \text{where}\ \ D_{n}=\frac{1}{2\pi i}\int_{C^{-}}D(z)z^{-n-1}\text{d}z.$
(3.16)

With this Laurent series, we can easily perform the sum-split,
$D(z)=D^{+}(z)+D^{-}(z)$,

$\displaystyle D^{+}(z)=\sum_{n=0}^{\infty}D_{n}z^{n},\ \
D^{-}(z)=\sum_{n=-\infty}^{-1}D_{n}z^{n}.$ (3.17)

This results in the final Wiener–Hopf equation,

$\displaystyle K^{+}(z)A^{+}(z)$
$\displaystyle-\frac{F_{\text{pole}}^{+}(z)}{K^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})}-D^{+}(z)$
(3.18) $\displaystyle=\ $ $\displaystyle
F_{\text{pole}}^{+}(z)\left[\frac{1}{K^{-}(z)}-\frac{1}{K^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})}\right]+D^{-}(z)+\frac{F^{-}(z)}{K^{-}(z)},$
(3.19)

where the left-hand (resp. right-hand) side is analytic on $\Omega^{+}$ (resp.
$\Omega^{-}$). The two sides of this Wiener–Hopf equation are used to
construct a function that is entire on the $z$ complex plane,

$\displaystyle\Psi(z)=\begin{cases}\eqref{HW-A-WHE-final-
LHS}&z\in\Omega^{+},\\\ \eqref{HW-A-WHE-final-RHS}&z\in\Omega^{-},\\\
\eqref{HW-A-WHE-final-LHS}=\eqref{HW-A-WHE-final-
RHS}&z\in\Omega^{+}\cap\Omega^{-}.\\\ \end{cases}$ (3.20)

Each term on the right-hand side is $O\left(\frac{1}{z}\right)$ as
$|z|\rightarrow\infty$, which implies that $\Psi$ is bounded and tends to zero
at infinity as well as entire. Therefore, Liouville’s theorem ensures that
both (3.18) and (3.19) are equivalently zero, and then we obtain the solution
for $A^{+}(z)$,

$\displaystyle
A^{+}(z)=\frac{F_{\text{pole}}^{+}(z)}{K^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})K^{+}(z)}+\frac{D^{+}(z)}{K^{+}(z)}.$
(3.21)

To get the scattering coefficients, we must apply the inverse Z-transform to
(3.21),

$\displaystyle A_{m}=$ $\displaystyle\
\frac{e^{iks\cos(\theta_{\text{I}}-\alpha)}}{2\pi
iK^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})}\oint_{C^{-}}\frac{z^{-m-1}}{K^{+}(z)(z-e^{iks\cos(\theta_{\text{I}}-\alpha)})}\text{d}z$
(3.22) $\displaystyle+\frac{1}{2\pi
i}\oint_{C^{-}}\frac{D^{+}(z)}{K^{+}(z)}z^{-m-1}\text{d}z.$

The integration contour $C^{-}$ takes the path that runs around the unit
circle with $z=0$ being the only potential pole inside. The left integral in
(3.22) is the solution for the top array in isolation $A^{(0)}_{m}$, which is
equivalent to the semi-infinite array solution $A^{\text{(s-inf)}}_{m}$ (see
(2.20) and (2.21)), that is

$\displaystyle\frac{e^{iks\cos(\theta_{\text{I}}-\alpha)}}{2\pi
iK^{-}(e^{iks\cos(\theta_{\text{I}}-\alpha)})}\oint_{C^{-}}\frac{z^{-m-1}}{K^{+}(z)(z-e^{iks\cos(\theta_{\text{I}}-\alpha)})}\text{d}z=A^{(0)}_{m}.$
(3.23)

To evaluate the right integral in (3.22), we should recall the Laurent series,
(2.23) and (3.17),

$\displaystyle\frac{1}{2\pi
i}\oint_{C^{-}}\frac{D^{+}(z)}{K^{+}(z)}z^{-m-1}\text{d}z=\sum_{p=0}^{\infty}\sum_{n=0}^{\infty}\frac{\lambda_{p}D_{n}}{2\pi
i}\oint_{C^{-}}z^{p+n-m-1}\text{d}z.$ (3.24)

The value of this integral is non-zero only when $p+n-m=0$. This implies that
the only non-zero value of the $p$-summation is $p=m-n\geq 0$ which implies
that,

$\displaystyle\frac{1}{2\pi
i}\oint_{C^{-}}\frac{D^{+}(z)}{K^{+}(z)}z^{-m-1}\text{d}z=\sum_{n=0}^{m}\lambda_{m-n}D_{n}.$
(3.25)

The coefficients $D_{n}$ are given by (3.16),

$\displaystyle D_{n}$
$\displaystyle=-\sum_{q=1}^{\infty}\sum_{p=0}^{\infty}\sum_{l=0}^{\infty}\lambda_{p}B_{-q}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(q,l)\right)\frac{1}{2\pi
i}\int_{C^{-}}z^{l-n-p-1}\text{d}z,$ (3.26)

where the integral is non-zero only when $l-n-p=0$ which implies that,

$\displaystyle D_{n}$
$\displaystyle=-\sum_{q=1}^{\infty}\sum_{p=0}^{\infty}\lambda_{p}B_{-q}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(q,p+n)\right).$
(3.27)

This means that the scattering coefficients $A_{m}$ for all $m\geq 0$ are
written as,

$\displaystyle A_{m}=A^{(0)}_{m}$
$\displaystyle-\sum_{q=1}^{\infty}\sum_{p=0}^{\infty}\sum_{n=0}^{m}\lambda_{m-n}\lambda_{p}B_{-q}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(q,p+n)\right),\quad
m\geq 0,$ (3.28)

in terms of the coefficients $\lambda_{m}$ given by (2.24) and the bottom
array coefficients $B_{-m}$. Note that setting $m<0$ in (3.28) automatically
gives $A_{m}=0$ as required.

###  3.2 Discrete Wiener–Hopf technique on bottom array scattering
coefficients

The other system of equations (3.6) can also be solved using the discrete
Wiener–Hopf technique for semi-infinite arrays. As before, we extend (3.6)
using some unknown coefficients $F_{m}$ and recalling $B_{m}=0$ for $m\geq 0$,

$\displaystyle B_{m}H^{(1)}_{0}(ka)$
$\displaystyle+\sum_{\begin{subarray}{c}n=-\infty\\\ n\neq
m\end{subarray}}^{\infty}\\!B_{n}H^{(1)}_{0}\\!(ks|m-n|)$ (3.29)
$\displaystyle=\begin{dcases}F_{m},&m\geq 0,\\\
-\sum_{n=0}^{\infty}\\!A_{n}H^{(1)}_{0}\\!\left(ks\Lambda_{\alpha}(-m,n)\right)\\!-e^{iksm\cos(\theta_{\text{I}}+\alpha)},&m<0.\end{dcases}$

Contrary to (3.7), the $B_{m}$ scattering coefficients are the unknowns to
find and the $A_{m}$ coefficients are assumed to be known. We also use a
slight variation of the Z-transform to (3.8) in this problem,

$\displaystyle G(z)=\sum_{m=-\infty}^{\infty}G_{m}z^{-m-1},\quad
G_{m}=\frac{1}{2\pi i}\oint_{C^{-}}G(z)z^{m}\text{d}z.$ (3.30)

This new Z-transform is applied to (3.29) to produce the Wiener–Hopf equation,

$\displaystyle K(z)B^{+}(z)=$ $\displaystyle
F_{\text{pole}}^{+}(z)+F_{A}^{+}(z)+F^{-}(z),$ (3.31)

where $B^{+}(z)$ is the transform of the bottom array scattering coefficients,

$\displaystyle
B^{+}(z)=\sum_{m=-\infty}^{\infty}B_{m}z^{-m-1}=\sum_{m=-\infty}^{-1}B_{m}z^{-m-1}=\sum_{m=0}^{\infty}B_{-m-1}z^{m}.$
(3.32)

This problem reuses the $\pm$ superscript notation and the definitions of the
analytic regions $\Omega^{\pm}$, the Wiener–Hopf kernel $K(z)$ (2.15) and the
factorisations $K^{\pm}(z)$. The three forcing terms on the right-hand side of
the Wiener–Hopf equation are defined as,

$\displaystyle F_{\text{pole}}^{+}(z)$
$\displaystyle=-\sum_{m=-\infty}^{-1}e^{iksm\cos(\theta_{\text{I}}+\alpha)}z^{-m-1}=\frac{1}{z-e^{iks\cos(\theta_{\text{I}}+\alpha)}},$
(3.33) $\displaystyle F_{A}^{+}(z)$
$\displaystyle=-\sum_{m=1}^{\infty}\sum_{n=0}^{\infty}A_{n}z^{m-1}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(m,n)\right),$
(3.34) $\displaystyle F^{-}(z)$
$\displaystyle=\sum_{m=0}^{\infty}F_{m}z^{-m-1}.$ (3.35)

We proceed with the discrete Wiener–Hopf technique in the same way as in
Section 3.1 and apply the new inverse Z-transform to determine $B_{m}$. In the
end, the bottom array scattering coefficients are given by,

$\displaystyle B_{m}=B^{(0)}_{m}$
$\displaystyle-\sum_{q=0}^{\infty}\sum_{p=0}^{\infty}\sum_{n=1}^{-m}\lambda_{-m-n}\lambda_{p}A_{q}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(q,p+n)\right),\quad
m<0,$ (3.36)

in terms of the coefficients $\lambda_{m}$ (2.24) and the top array
coefficients $A_{m}$. Note that setting $m\geq 0$ in (3.36) automatically
gives $B_{m}=0$ as required.

###  3.3 Point scatterer wedge iterative scheme

The two formulae for the scattering coefficients given by (3.28) and (3.36)
are still coupled and are difficult to solve exactly. However, we can use them
to form an iterative scheme with the isolated solution $A^{(0)}_{m}$ and
$B^{(0)}_{-m}$ as an initial guess

$\displaystyle A^{(0)}_{m}$
$\displaystyle=-\frac{e^{-iksm\cos(\theta_{\text{I}}-\alpha)}}{K^{+}(e^{-iks\cos(\theta_{\text{I}}-\alpha)})}\sum_{n=0}^{m}\lambda_{n}e^{iksn\cos(\theta_{\text{I}}-\alpha)},$
(3.37) $\displaystyle B^{(0)}_{-m}$
$\displaystyle=-\frac{e^{-iksm\cos(\theta_{\text{I}}+\alpha)}}{K^{+}(e^{-iks\cos(\theta_{\text{I}}+\alpha)})}\sum_{n=0}^{m-1}\lambda_{n}e^{iksn\cos(\theta_{\text{I}}+\alpha)}.$

Let us denote $A^{(j)}_{m}$ and $B^{(j)}_{-m}$ as the $j^{\text{th}}$
approximations to the scattering coefficients $A_{m}$ and $B_{-m}$
respectively. For the first iteration, we substitute the initial guess
$B^{(0)}_{-m}$ into (3.28) to get the first iteration of the top array
scattering coefficients $A^{(1)}_{m}$,

$\displaystyle
A^{(1)}_{m}=A^{(0)}_{m}-\sum_{q=1}^{\infty}\sum_{p=0}^{\infty}\sum_{n=0}^{m}\\!\lambda_{m-n}\lambda_{p}B^{(0)}_{-q}H^{(1)}_{0}\\!\left(ks\Lambda_{\alpha}(q,p+n)\right).$

Next, we substitute $A^{(1)}_{m}$ into (3.36) to get the first iteration of
the bottom array coefficients $B^{(1)}_{-m}$,

$\displaystyle
B^{(1)}_{-m}=B^{(0)}_{-m}-\sum_{q=0}^{\infty}\sum_{p=0}^{\infty}\sum_{n=1}^{m}\\!\lambda_{m-n}\lambda_{p}A^{(1)}_{q}H^{(1)}_{0}\\!\left(ks\Lambda_{\alpha}(q,p+n)\right).$

We can repeat this process as many times as needed to get the $j^{\text{th}}$
approximations, $A^{(j)}_{m}$ and $B^{(j)}_{-m}$ with the following general
equations,

$\displaystyle A^{(j)}_{m}$
$\displaystyle=A^{(0)}_{m}-\sum_{q=1}^{\infty}\sum_{p=0}^{\infty}\sum_{n=0}^{m}\\!\lambda_{m-n}\lambda_{p}B^{(j-1)}_{-q}H^{(1)}_{0}\\!\left(ks\Lambda_{\alpha}(q,p+n)\right),$
$\displaystyle B^{(j)}_{-m}$
$\displaystyle=B^{(0)}_{-m}-\sum_{q=0}^{\infty}\sum_{p=0}^{\infty}\sum_{n=1}^{m}\\!\lambda_{m-n}\lambda_{p}A^{(j)}_{q}H^{(1)}_{0}\\!\left(ks\Lambda_{\alpha}(q,p+n)\right).$

Consider the two infinite vectors $\mathbf{A}^{(j)}$ and $\mathbf{B}^{(j)}$
which include the $j^{\textrm{th}}$ iteration of the scattering coefficients
from $A^{(j)}_{0}$ to $A^{(j)}_{\infty}$ and $B^{(j)}_{-1}$ to
$B^{(j)}_{-\infty}$ respectively. We can write the iterative scheme as the
following infinite matrix system of equations.

$\displaystyle\mathbf{A}^{(j)}$
$\displaystyle=\mathbf{A}^{(0)}-\mathcal{M}^{(B)}\mathbf{B}^{(j-1)},$
$\displaystyle\mathbf{B}^{(j)}$
$\displaystyle=\mathbf{B}^{(0)}-\mathcal{M}^{(A)}\mathbf{A}^{(j)},$ (3.38)

where the matrices $\mathcal{M}^{(A)}$ and $\mathcal{M}^{(B)}$ have the
following entries,

$\displaystyle\mathcal{M}^{(A)}_{mq}$
$\displaystyle=\sum_{p=0}^{\infty}\sum_{n=1}^{m}\lambda_{m-n}\lambda_{p}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(q,p+n)\right),\
\ (m>0,\ q\geq 0),$ $\displaystyle\mathcal{M}^{(B)}_{mq}$
$\displaystyle=\sum_{p=0}^{\infty}\sum_{n=0}^{m}\lambda_{m-n}\lambda_{p}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(q,p+n)\right),\
\ (m\geq 0,\ q>0).$

In practice, we need to truncate the infinite sums at a sufficiently high
point. For example, say that we truncate both sums at $M$, then we need to
evaluate $M+1$ $\lambda_{n}$ coefficients (i.e. $\lambda_{0},\ \lambda_{1},\
...\ \lambda_{M}$), for which we use (A.17). This gives $\mathcal{M}^{(A)}$ as
a $M\times(M+1)$ matrix and $\mathcal{M}^{(B)}$ as a $(M+1)\times M$ matrix.
Via (3.37), that same list of $\lambda_{n}$ also gives us $M+1$ $A^{(0)}_{m}$
and $M$ $B^{(0)}_{-m}$ coefficients (i.e. $A^{(0)}_{0}$, … $A^{(0)}_{M}$ and
$B^{(0)}_{-1}$, … $B^{(0)}_{-M}$) as an initial guess. With the iteration
matrices and the initial guess defined, we then apply (3.3) $j$ times to get
$M+1$ $A^{(j)}_{m}$ and $M$ $B^{(j)}_{-m}$ coefficients at each stage of the
iterative scheme.

###  3.4 Convergence of the iterative scheme

There are several different parameters that can affect the conditions for
convergence of the iterative scheme: $\alpha$, $ks$, $ka$ and the sum
truncation numbers. To assess the convergence conditions, we define the
$j^{\text{th}}$ error vector of the scattering coefficients by
$\mathbf{E}^{(j)}_{A}=\mathbf{A}^{(j)}-\mathbf{A}$ and
$\mathbf{E}^{(j)}_{B}=\mathbf{B}^{(j)}-\mathbf{B}$ which satisfies the
following matrix system,

$\displaystyle\mathbf{E}^{(j)}_{A}$
$\displaystyle=-\mathcal{M}^{(B)}\mathbf{E}^{(j-1)}_{B},$
$\displaystyle\mathbf{E}^{(j)}_{B}$
$\displaystyle=-\mathcal{M}^{(A)}\mathbf{E}^{(j)}_{A}.$ (3.39)

We can substitute each equation of (3.4) into the other to obtain the
following,

$\displaystyle\mathbf{E}^{(j)}_{A}$
$\displaystyle=\mathcal{M}^{(B)}\mathcal{M}^{(A)}\mathbf{E}^{(j-1)}_{A},$
$\displaystyle\mathbf{E}^{(j)}_{B}$
$\displaystyle=\mathcal{M}^{(A)}\mathcal{M}^{(B)}\mathbf{E}^{(j-1)}_{B}.$
(3.40)

It is well known that the error vectors $\mathbf{E}^{(j)}_{A}$ and
$\mathbf{E}^{(j)}_{B}$ will tend to zero if and only if the spectral radius of
either iteration matrices $\mathcal{M}^{(A)}\mathcal{M}^{(B)}$ and
$\mathcal{M}^{(B)}\mathcal{M}^{(A)}$ is less than 1, which implies that the
iterative scheme will converge.

We plot the spectral radius of $\mathcal{M}^{(A)}\mathcal{M}^{(B)}$ (and call
it $\rho$) in Figure 4 with respect to the wedge angle $\alpha$. It is
important to note that $\rho$ is symmetric about $\alpha=\frac{\pi}{2}$, and
although we have plotted $\alpha\in(0,\ \pi)$, the true range of interest here
should be $\alpha\in(\sin^{-1}(a/s),\ \pi-\sin^{-1}(a/s))$ which satisfies the
condition $\sin(\alpha)>a/s$ (recall this is to prevent overlapping cylinders
when $\alpha\approx 0$ or $\pi$). Figure 4 shows that all the considered cases
in this $\alpha$ range will therefore imply that the iterative scheme will
converge to a solution.

Figure 4: Here is the spectral radius $\rho$ of the matrix
$\mathcal{M}^{(A)}\mathcal{M}^{(B)}$ with respect to the angle $\alpha$ and
summations truncated at $M=1000$. Each plot on the left side has a different
value for the wavenumber $k$ while preserving $ka\ll 2\pi$. Similarly, the
right side varies the cylinder radius $a$ within the constraint $a<s/2$. Note
that the horizontal dashed line indicates the critical spectral radius above
which the iterative scheme will stop converging. We also mark the edges of the
condition $\sin(\alpha)>a/s$ by vertical dashed lines (left) or crosses
(right).

In producing these graphs, we have observed that $\rho$ converges very quickly
with respect to the sum truncation number $M$, so we have deemed it sufficient
to set $M=1000$. Figure 4 (left) illustrates the effect of the wavenumber $k$
on the spectral radius. It shows that the spectral radius is smooth with
respect to $\alpha$ when $ks<\pi$, but becomes increasingly jagged as $k$ gets
larger (this might be alleviated by increasing $M$). The exception to this is
when $ks$ is a multiple of $\pi$, which is a special case that we want to
study in a future article.

Figure 4 (right) displays how the cylinder radius $a$ affects the convergence.
It clearly shows that the iterative scheme will converge faster as $a$ becomes
smaller. This trend also continues for larger wavenumbers, however, $a\ll s/2$
will need to be more strictly satisfied to keep $\rho$ below the critical
value.

For higher values of $k$ and $s$, the plots in Figure 4 will appear to be less
smooth with respect to $\alpha$. This is due to a change in dominance between
the eigenvalues that produce the spectral radius. Apart from some rare
exceptions, we can reduce the overall non-smoothness by increasing the value
of the sum truncation number $M$, but it should be noted that the computation
time is of the order $O(M^{3})$.

Repeating the same analysis for the other iteration matrix
$\mathcal{M}^{(B)}\mathcal{M}^{(A)}$ yields almost identical results (the
observed difference is so small that rounding error is the most likely cause).
Therefore, we are satisfied that in most cases, the iterative scheme will
converge provided that our initial assumptions from Foldy’s approximation
(namely $a\ll s/2$ and $ka\ll 2\pi$) are reasonably satisfied. Additionally,
the iterative scheme will perform at its best when $ks<\pi$ as well.

###  3.5 Numerical results of the wave field

Now that we have assessed the spectral radius of the iteration matrices, we
have a good idea on what parameters will lead to a successful convergence.
However, we have not considered the special case of resonance. Referred to as
Wood’s anomalies after being observed experimentally in Wood (1902), this
phenomenon occurs when the scattered waves of each unit cell are synchronised
and pieced together to form constructive interference as a consequence. For
infinite and semi-infinite arrays, the total scattered field usually takes the
form of a single plane wave propagating parallel to the array. For the wedge,
there are four conditions for resonance (two per array) given by,

$\displaystyle\frac{ks}{2\pi}\left(1\mp\cos(\psi-\alpha)\right)\in{\mathbb{Z}},\
\ \text{and}\ \
\frac{ks}{2\pi}\left(1\mp\cos(\psi+\alpha)\right)\in{\mathbb{Z}},$ (3.41)

where $\psi$ could be the propagation angle of the incident wave or any of the
reflected or transmitted wave angles, as determined from geometrical optics.
This is a very interesting special case to study and we shall do so in a
future article. In the following, we shall choose test cases that are
physically viable and non-resonance cases and then, once we are satisfied with
the approximations $A^{(j)}_{m}$ and $B^{(j)}_{-m}$, we substitute them into
(3.2) to get the physical wave fields $\Phi_{\text{S}}$ and $\Phi$.

For two of these test cases, we perform fifty iterations and, in Figure 5,
plot the infinity norm of the difference between the $j^{\textrm{th}}$
iteration and the final iteration of the scattering coefficients. In other
words, we plot the quantities
$||\mathbf{A}^{(j)}-\mathbf{A}^{(50)}||_{\infty}$ and
$||\mathbf{B}^{(j)}-\mathbf{B}^{(50)}||_{\infty}$ in red and blue
respectively. Note that the vectors here include coefficients up to and
including $A^{(j)}_{100}$ and $B^{(j)}_{-100}$. The magenta line is given by
$\rho^{j}$ where $\rho$ is the spectral radius, which shows that all these
errors have the order $O(\rho^{j})$. For these test cases, this figure shows
that the iteration scheme is becoming accurate to machine precision by the
$25^{\textrm{th}}$ iteration.

Figure 5: We plot the infinity norm error of the scattering coefficients.
Here, the red (resp. blue) line displays the error quantity
$||\mathbf{A}^{(j)}-\mathbf{A}^{(50)}||_{\infty}$ (resp.
$||\mathbf{B}^{(j)}-\mathbf{B}^{(50)}||_{\infty}$) with respect to iteration
stage $j$, where the vectors contain the coefficients up to and including
$A^{(j)}_{100}$ (resp. $B^{(j)}_{-100}$) The magenta line is the order of the
convergence given by the spectral radius $\rho$. The summations are truncated
at $M=1000$, the wedge is defined by $\alpha=\frac{5\pi}{6}$, $s=0.1$ and
$a=0.01$, and the left (resp. right) plot has an incident wave with the
parameters $k=5\pi$ and $\theta_{\textrm{I}}=0$ (resp. $k=15\pi$ and
$\theta_{\textrm{I}}=\frac{\pi}{2}$).

Figure 6 features the scattered wave field’s real part of two different test
cases using the iterative scheme and compares them to the equivalent finite-
element simulations determined through COMSOL. In these figures, we iterated
the scattering coefficients twenty-five times and truncated the summations at
$M=1000$. This corresponds to considering the interaction of 1000 cylinders on
both the top and bottom face as well as the one at the wedge edge. The COMSOL
simulations are limited by the number of installed scatterers, as more
scatterers requires a larger numerical mesh to find an accurate solution. In
these simulations, we created the point scatterer wedge from a total of 61
cylinders and used perfectly-matched layers to prevent incoming radiation.
Between the left and right sides of Figure 6, we can clearly see that the
agreement is excellent. We have found that in the numerous other test cases we
compared, we similarly have excellent agreement.

Figure 6: Here, we plot the real part of the scattered wave field where the
point scatterer wedge is defined by $\alpha=\frac{5\pi}{6}$, $s=0.1$ and
$a=0.01$. The left column is using the iterative scheme (twenty-five
iterations were performed and summations were truncated at $M=1000$) and the
right column plots are COMSOL simulations. The top row has an incident wave
with the parameters $k=5\pi$ and $\theta_{\textrm{I}}=0$. The bottom row has
the parameters $k=15\pi$ and $\theta_{\textrm{I}}=\frac{\pi}{2}$.

To see the evolution of the wave field as the iterative scheme progresses, we
plot the difference in the wave field $\Phi$ between the iterative scheme and
the semi-infinite array approximation. Specifically, Figure 7 (left) displays
the real part difference between Figure 6 (top left) and Figure 3 (left) which
shows this difference closely resembling a cylindrical wave emanating from the
region of the wedge edge. Figure 7 (right) takes
$\textup{Re}\left\\{\Phi(3s/2,\theta)\right\\}$ at several $\theta$ values and
illustrates the error between the first ten iterations and the final
$25^{\textrm{th}}$ iteration. It shows that, at most points, the convergence
is oscillatory. Further plots (not included here) of the absolute error give
the convergence order as $O(\rho^{j})$, the same behaviour as shown in Figure
5.

Figure 7: The left figure plots the real part difference in the wave field
between the $25^{\textrm{th}}$ iteration and the semi-infinite array
approximation. The right figure shows the real part error progression of
several wave field points. Here, the point scatterer wedge is defined by
$\alpha=\frac{5\pi}{6}$, $s=0.1$ and $a=0.01$ and the incident wave has
parameters $k=5\pi$ and $\theta_{\textrm{I}}=0$.

##  4 Conclusions and further work

To summarise, we have decomposed the two semi-infinite arrays in the point
scatterer wedge into two separate semi-infinite array problems and solved both
using the discrete Wiener–Hopf technique. Then, we took this solution and
constructed an iterative scheme to approximate the scattering coefficients.
Lastly, we assessed the convergence of the iterative scheme via the spectral
radius of the iteration matrices and did a side-by-side comparison of some
test cases with COMSOL simulations.

We concluded from Figure 4 that the iterative scheme works very well provided
that $\sin(\alpha)>\frac{a}{s}$ is satisfied, we are away from resonance and
the assumptions associated with Foldy’s approximation ($a\ll\frac{s}{2}$ and
$ka\ll 2\pi$) are sufficiently satisfied (in particular, the best convergence
occurs when $ks<\pi$). Figure 5 showed that the iterative scheme has a
convergence order of $O(\rho^{j})$ with spectral radius $\rho$ and number of
iterations $j$. We also concluded from Figure 6 that the iterative scheme
agrees closely with the COMSOL comparison and Figure 7 illustrated how the
wave field converges to a solution with respect to the iterations. We see that
the first iteration is adding the coupling effect between the two semi-
infinite arrays, while each subsequent iteration corrects this coupling effect
towards convergence. We also see that the coupling effect is likely to be
strongest at the wedge edge and this convergence is most likely exponentially
decaying and oscillatory.

Although we mentioned earlier that we deliberately chose test cases that do
not result in resonance, this does not mean that we are not interested in the
subject. Taking a naive approach, one can say that there are four conditions
for resonance which are given by,

$\displaystyle\frac{ks}{2\pi}(1\mp\cos(\alpha-\theta_{\textrm{I}}))\in{\mathbb{Z}},\quad\frac{ks}{2\pi}(1\mp\cos(\alpha+\theta_{\textrm{I}}))\in{\mathbb{Z}},$
(4.1)

two for the top array and two for the bottom array, respectively. However,
this does not take in account the interactions between the two arrays. For
example, could any of the scattered waves (or even resonance waves) from one
array cause resonance in the other one? If so, then there will be extra
conditions where $\theta_{\textrm{I}}$ in (4.1) is appropriately replaced for
each of these scattered waves. This means that this coupling effect can become
convoluted very quickly without adequate ray-tracing techniques to keep track
of these scattered waves and their reflections/transmissions. It is very
unlikely that the coupling effect influences the nature of the resonance waves
themselves, i.e. they will still be in the form of plane waves, but this is
unproven.

For the iterative scheme, there are some avenues which can improve the
accuracy and speed of numerical computations. These avenues will be especially
useful in parameter regimes with slow convergence. One particular avenue is
the method of tail-end asymptotics, which could be used to reduce the error of
truncations in infinite summations. For example, see Appendix A.1 where we
applied it to the Wiener-Hopf kernel $K(z)$. This method is also beneficial
for calculating both the spectral radius and the scattering coefficients.
However, in the latter, it requires prior knowledge of the asymptotics of the
scattering coefficients $A_{m}$ and $B_{-m}$ as $m\rightarrow\infty$ which can
be done with the same ray-tracing techniques mentioned earlier for tracking
resonance cases.

Alternatively, this problem can be approached by separating it into distinct
symmetric (labelled ‘even’) and anti-symmetric (labelled ‘odd’) sub-problems.
In these sub-problems, the incident field is given by,

$\displaystyle\Phi_{\text{I}}^{\textrm{(even)}}(r,\theta)=\frac{1}{2}\left(e^{-ikr\cos(\theta-\theta_{\text{I}})}+e^{-ikr\cos(\theta+\theta_{\text{I}})}\right),$
(4.2)
$\displaystyle\Phi_{\text{I}}^{\textrm{(odd)}}(r,\theta)=\frac{1}{2}\left(e^{-ikr\cos(\theta-\theta_{\text{I}})}-e^{-ikr\cos(\theta+\theta_{\text{I}})}\right),$
(4.3)

and the scattering coefficients satisfy the relations
$B_{-n}^{\textrm{(even)}}=A_{n}^{\textrm{(even)}}$ and
$B_{-n}^{\textrm{(odd)}}=-A_{n}^{\textrm{(odd)}}$. Each sub-problem leads to a
single system of equations to solve:

$\displaystyle A^{\textrm{(even)}}_{m}H^{(1)}_{0}(ka)$
$\displaystyle+\sum_{\begin{subarray}{c}n=0\\\ n\neq
m\end{subarray}}^{\infty}\left[A^{\textrm{(even)}}_{n}H^{(1)}_{0}(ks|m-n|)\right]$
$\displaystyle+\sum_{n=1}^{\infty}\left[A^{\textrm{(even)}}_{n}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(m,n)\right)\right]+\Phi_{\text{I}}^{\textrm{(even)}}(sm,\alpha)=0,$
$\displaystyle A^{\textrm{(odd)}}_{m}H^{(1)}_{0}(ka)$
$\displaystyle+\sum_{\begin{subarray}{c}n=0\\\ n\neq
m\end{subarray}}^{\infty}\left[A^{\textrm{(odd)}}_{n}H^{(1)}_{0}(ks|m-n|)\right]$
$\displaystyle-\sum_{n=1}^{\infty}\left[A^{\textrm{(odd)}}_{n}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(m,n)\right)\right]+\Phi_{\text{I}}^{\textrm{(odd)}}(sm,\alpha)=0,$
(4.4)

which are similar to (3.4) and (3.6), but this time the two systems are
decoupled. This means that if each system can be solved, then no iteration
procedure would be needed. However, each system would only be solvable using
the discrete Wiener-Hopf technique provided that we can explicitly write the
$Z$-transform of
$\sum_{n=1}^{\infty}\left[A_{n}H^{(1)}_{0}\left(ks\Lambda_{\alpha}(m,n)\right)\right]$
in the form of $F(z)A(z)$ where $F(z)$ will become an extra term in the
discrete Wiener-Hopf kernel. Since $\Lambda_{\alpha}(m,n)$ does depend on
$\alpha$ and not solely on $m-n$, we cannot use the discrete convolution
theorem to do this directly. Hence, at this stage, it is unknown what the
explicit formula for $F(z)$ would be. If it was known, then the rest of the
procedure should remain exactly the same as the semi-infinite array problem.
In the end, the technique would give an exact solution to each sub-problem and
then their sum will be the exact solution to the full problem.

In this article, we considered the simplest possible unit cell in order to
highlight the effects of edges interacting, but certainly many more
configurations are possible. Much more complicated models of array scattering
with Helmholtz resonators can likewise be considered within the framework of
Foldy’s approximation Brandão and Schnitzer (2020). We also assumed that the
cylinder scatterers were isotropic and hence, restricted the study to sound-
soft (or Dirichlet) boundary conditions. However, there are other boundary
types that we can consider such as sound-hard (Neumann), impedance or
penetrable scatterers for example. In these cases, the isotropic assumption
made for Foldy’s approximation will not be sufficient. Instead, we will need
to at least include double pole terms in the multi-pole expansion which are
represented by the Hankel functions $H^{(1)}_{1}(z)$ (see section 8.3.3 in
Martin (2006)). This means that we would have to replace the general solution
(2.2) with,

$\displaystyle\Phi_{\text{S}}(\mathbf{r})$
$\displaystyle=\sum_{n}\left[A_{n}H^{(1)}_{0}(kr_{n})+(a_{n}\cos(\theta_{n})+b_{n}\sin(\theta_{n}))H^{(1)}_{1}(kr_{n})\right],$
(4.5)

where $r_{n}=|\mathbf{r}-\mathbf{R}_{n}|$ and
$\theta_{n}=\arg(\mathbf{r}-\mathbf{R}_{n})$. In this new form, each local
scattered wave features the original monopole term as well as an additional
dipole term. However, for different boundary conditions, we anticipate that
additional difficulties will arise. This is due to the fact that infinite and
semi-infinite arrays with Neumann boundary conditions can support Rayleigh-
Bloch surface waves, as discussed in depth in Linton et al. (2007).

Apart from Foldy’s approximation, another interesting mathematical approach at
low frequency is that of asymptotic homogenisation. In particular, for
straight or curved arrays of scatterers, the result of such approach is a
continuous boundary with effective jump conditions on both the field and its
normal derivative Chapman et al. (2015); Hewett and Hewitt (2016); Touboul et
al. (2020a, b). Such boundary is referred to as a meta-interface in the
literature. In our case, the meta-interface will be that of a wedge and will
hence include a geometric singularity at the corner. It will be very
interesting indeed to develop mathematical and/or numerical methods to tackle
this ‘homogenised’ point scatterer wedge problem and compare the results to
our existing approach.

Although we designed this iterative scheme with a specific arrangement of
scatterers in mind, it is in principle possible to use it with two semi-
infinite arrays that do not share a common edge and/or have different spacing
parameters ($s_{1}$ and $s_{2}$ say). Let us fix one of the arrays to be the
traditional semi-infinite array supported on the positive $x$-axis (with
$s=s_{1}$). Then the form of the function $\Lambda_{\alpha}(m,n)$ will be
determined from the relative positioning of the scatterers between the two
arrays and will depend on the angle ($\alpha$), the spacing ($s_{2}$) and the
origin of the second array. The only additional necessary condition here for
the iterative method to work, is that no two scatterers overlap or touch in
any way (as we did for the wedge in this article). Provided that this
condition is met, the methodology of this article should be unchanged. Whether
or not the iterative scheme converges will be determined by the resulting
iterative matrices and their spectral radii.

## References

  * Bennetts et al. [2017] L. G. Bennetts, M. A. Peter, and F. Montiel.  Localisation of Rayleigh-Bloch waves and damping of resonant loads on arrays of vertical cylinders.  _J. Fluid Mech._ , 813:508–527, 2017.
  * Brandão and Schnitzer [2020] R. Brandão and O. Schnitzer.  Asymptotic modeling of Helmholtz resonators including thermoviscous effects.  _Wave Motion_ , 97(102583):1–25, 2020.
  * Bruno and Fernandez-Lado [2017] O. P. Bruno and A. G. Fernandez-Lado.  Rapidly convergent quasi-periodic Green functions for scattering by arrays of cylinders—including Wood anomalies.  _Proc. Roy. Soc. A_ , 473(2199):20160802, 2017\.
  * Camacho et al. [2019] M. Camacho, A. P. Hibbins, F. Capolino, and M. Albani.  Diffraction by a truncated planar array of dipoles: A Wiener–Hopf approach.  _Wave Motion_ , 89:28–42, 2019.
  * Chapman et al. [2015] S. J. Chapman, D. P. Hewett, and L. N. Trefethen.  Mathematics of the Faraday cage.  _SIAM Rev._ , 57(3):398–417, 2015.
  * Craster and Guenneau [2013] R. V. Craster and S. Guenneau.  _Acoustic Metamaterials_.  Springer, 2013.
  * Foldy [1945] L. L. Foldy.  The multiple scattering of waves. I. General theory of isotropic scattering by randomly distributed scatterers.  _Phys. Rev._ , 67(3-4):107–119, 1945.
  * Gradshteyn and Ryzhik [2014] I. S. Gradshteyn and I. M. Ryzhik.  _Table of Integrals, Series, and Products_.  Academic Press, 8th edition, 2014.
  * Hewett and Hewitt [2016] D. P. Hewett and I. J. Hewitt.  Homogenized boundary conditions and resonance effects in Faraday cages.  _Proc. Roy. Soc. A_ , 472(2189), 2016.
  * Hills and Karp [1965] N. L. Hills and S. N. Karp.  Semi-Infinite Diffraction Gratings-I.  _Comm. Pure Appl. Math_ , XVIII:203–233, 1965.
  * Kamotskiĭ and Nazarov [1999] I. V. Kamotskiĭ and S. A. Nazarov.  Wood’s anomalies and surface waves in the problem of scattering by a periodic boundary. I.  _Mat. Sb._ , 190(1):109–138, 1999.
  * Korolkov et al. [2016] A. I. Korolkov, S. A. Nazarov, and A. V. Shanin.  Stabilizing solutions at thresholds of the continuous spectrum and anomalous transmission of waves.  _ZAMM Z. Angew. Math. Mech._ , 96(10):1245–1260, 2016.
  * Linton [1998] C. M. Linton.  The Green’s function for the two-dimensional Helmholtz equation in periodic domains.  _J. Eng. Math._ , 33(4):377–402, 1998.
  * Linton and Martin [2004] C. M. Linton and P. A. Martin.  Semi-Infinite Arrays of Isotropic Point Scatterers. A Unified Approach.  _SIAM J. Appl. Math_ , 64(3):1035–1056, 2004\.
  * Linton et al. [2007] C. M. Linton, R. Porter, and I. Thompson.  Scattering by a semi-infinite periodic array and the excitation of surface waves.  _SIAM J. Appl. Math_ , 67(5):1233–1258, 2007\.
  * Liu and Declercq [2015] J. Liu and N. F. Declercq.  Investigation of the origin of acoustic wood anomaly.  _J. Acoust. Soc. Am._ , 138(2):1168–1179, 2015\.
  * Lynott et al. [2019] G. M. Lynott, V. Andrew, I. D. Abrahams, M. J. Simon, W. J. Parnell, and R. C. Assier.  Acoustic scattering from a one-dimensional array; tail-end asymptotics for efficient evaluation of the quasi-periodic green’s function.  _Wave Motion_ , 89:232–244, 2019.
  * Martin [2006] P. A. Martin.  _Multiple Scattering: Interaction of Time-Harmonic Waves with N obstacles_.  Cambridge University Press, Cambridge, 2006.
  * Movchan et al. [2018] A. B. Movchan, N. V. Movchan, I. S. Jones, and D. J. Colquitt.  _Mathematical Modelling of Waves in Multi-Scale Structured Media_.  Chapman & Hall/CRC, 2018.
  * Nazarov [2008] S. A. Nazarov.  The Neumann problem in angular domains with periodic and parabolic perturbations of the boundary.  _Tr. Mosk. Mat. Obs._ , 69:182–241, 2008.
  * Nethercote et al. [2020] M. A. Nethercote, R. C. Assier, and I. D. Abrahams.  Analytical methods for perfect wedge diffraction: a review.  _Wave Motion_ , 93(102479), 2020.
  * Rayleigh [1907] L. Rayleigh.  On the Dynamical Theory of Gratings.  _Proc. Roy. Soc. A_ , 79:399–416, 1907.
  * Scarpetta and Sumbatyan [1996] E. Scarpetta and M. A. Sumbatyan.  Explicit analytical results for one-mode oblique penetration into a periodic array of screens.  _IMA J. Appl. Math._ , 56(2):109–120, 1996.
  * Shanin and Korolkov [2017] A. V. Shanin and A. I. Korolkov.  Diffraction of a mode close to its cut-off by a transversal screen in a planar waveguide.  _Wave Motion_ , 68:218–241, 2017.
  * Touboul et al. [2020a] M. Touboul, B. Lombard, and C. Bellis.  Time-domain simulation of wave propagation across resonant meta-interfaces.  _J. Comput. Phys_ , 414(109474), 2020a.
  * Touboul et al. [2020b] M. Touboul, K. Pham, A. Maurel, J. Marigo, B. Lombard, and C. Bellis.  Effective resonant model and simulations in the time-domain of wave scattering from a periodic row of highly-contrasted inclusions.  _J. Elast_ , 142:53–82, 2020b.
  * Wegert [2012] E. Wegert.  _Visual Complex Functions_.  Birkhauser Basel, 2012.
  * Wood [1902] R. W. Wood.  On a remarkable case of uneven distribution of light in a diffraction grating spectrum.  _Proc. Phys. Soc._ , 18(1):269–275, 1902.

##  Appendix A Asymptotics and numerical methods for the discrete Wiener–Hopf
kernel and factors

Infinite sums are required to be truncated repeatedly in this article. In this
appendix section, we shall discuss two different computational methods to
evaluate the discrete Wiener–Hopf kernel. The first is the method of tail-end
asymptotics to improve the convergence of the original formula. The second is
an alternative fast-convergent formula. Recall the original formula for the
discrete Wiener–Hopf kernel (2.15),

$\displaystyle
K(z)=H^{(1)}_{0}(ka)+\sum_{\ell=1}^{\infty}(z^{\ell}+z^{-\ell})H^{(1)}_{0}(ks\ell).$
(A.1)

In the limit as $\textup{Im}\left\\{k\right\\}\rightarrow 0$, this formula
only converges when $z$ is on the unit circle (not including the branch points
$z=e^{\pm iks}$) and even then, the convergence is conditional and very slow.
This is because the sum term has the following order as
$\ell\rightarrow\infty$,

$\displaystyle(z^{\ell}+z^{-\ell})H^{(1)}_{0}(ks\ell)=O\left((z^{\ell}+z^{-\ell})\frac{e^{iks\ell}}{\sqrt{\ell}}\right).$
(A.2)

Computing the kernel will naturally require us to truncate the sum at $L$ say.
This introduces an error term which we call the tail-end.

###  A.1 Method of tail-end asymptotics

Let us define the tail-end sum by,

$\displaystyle\text{Tail}(z,L)=\sum_{\ell=L}^{\infty}(z^{\ell}+z^{-\ell})H^{(1)}_{0}(ks\ell).$
(A.3)

If $L$ was sufficiently large, one could neglect the tail-end altogether, but
the slow convergence of (A.1) would make this impractical. Instead, the tail-
end sum is asymptotically approximated as $L\rightarrow\infty$ (provided that
$z\neq e^{\pm iks}$) using the same methodology as in Lynott et al. [2019] for
the quasi-periodic Green’s function. With two corrections, the result of this
method of tail-end asymptotics is given by,

$\displaystyle\text{Tail}(z,L)$ $\displaystyle\sim\frac{1-i}{\sqrt{\pi
ksL}}\left[T(ze^{iks},L)+T(e^{iks}/z,L)\right],$ (A.4)
$\displaystyle\text{where}\ T(z,L)$
$\displaystyle=\frac{z^{L}}{1-z}\bigg{[}1-\frac{1}{2L}\left(\frac{z}{1-z}+\frac{i}{4ks}\right)\bigg{]}.$

Even just one tail-end correction will drastically improve the convergence of
the Hankel summation formula. Specifically $N$ tail-end corrections will imply
that,

$\displaystyle K(z)=$ $\displaystyle\
H^{(1)}_{0}(ka)+\sum_{\ell=1}^{L-1}(z^{\ell}+z^{-\ell})H^{(1)}_{0}(ks\ell)+\text{Tail}(z,L)$
(A.5)
$\displaystyle+O\left(\left(\frac{z^{L}}{1-ze^{iks}}+\frac{z^{-L}}{1-e^{iks}/z}\right)\frac{e^{iksL}}{L^{\frac{1}{2}+N}}\right).$

To reinforce our statements, Figure 8 compares the various rates of
convergence of (A.5) at $z=1$ with a varying truncation point $L$ and
different numbers of tail-end corrections $N$. These plots clearly show that
no tail-end corrections will cause very slow convergence while even just one
tail-end correction is a significant improvement. Evidently, the method of
tail-end asymptotics is a powerful technique for slow-convergent infinite sums
like (A.1). It is also useful for improving the convergence of other infinite
sums in this article including the physical field (2.2) and the scattering
coefficients (2.8), (2.22), (3.28) and (3.36). However, these cases will
require us to determine the asymptotic behaviour of the $\lambda_{n}$
coefficients as $n\rightarrow\infty$.

###  A.2 Fast-convergent formula

Besides the original formula (A.1), we can reuse equations (8.522) and (8.524)
from Gradshteyn and Ryzhik [2014] to write the discrete Wiener–Hopf kernel in
a form that is faster to converge. Before we do this, it is convenient to map
$z$ onto a new complex plane using the following transformation $z=e^{it}$.
Under this mapping the $z$ unit circle is transformed onto the real line,
inside (resp. outside) the unit circle is mapped onto the upper (resp. lower)
half plane and (A.1) is rewritten as

$\displaystyle K(e^{it})=H^{(1)}_{0}(ka)+2\sum_{\ell=1}^{\infty}\cos(\ell
t)H^{(1)}_{0}(ks\ell).$ (A.6)

Applying equations (8.522) and (8.524) to (A.6) gives us the fast-convergent
formula,

$\displaystyle K(e^{it})=$ $\displaystyle\
H^{(1)}_{0}(ka)-1-\frac{2i}{\pi}\left(\gamma+\ln\left(\frac{ks}{4\pi}\right)\right)+\frac{2}{\sqrt{(ks)^{2}-t^{2}}}$
(A.7) $\displaystyle+\sum_{\begin{subarray}{c}\ell=-\infty\\\ \ell\neq
0\end{subarray}}^{\infty}\left(\frac{2}{\sqrt{(ks)^{2}-(2\pi\ell-t)^{2}}}+\frac{i}{\pi|\ell|}\right).$

where $\gamma=0.5772...$ is the Euler–Mascheroni constant. In this formula,
there are two sets of branch cuts starting at the points $t=\pm ks+2\pi\ell$.
The ones starting at $t=ks+2\pi\ell$ (resp. $t=-ks+2\pi\ell$) are orientated
to positive (resp. negative) imaginary infinity in the $t$ complex plane. All
of the square roots in (A.7) are defined to ensure that
$\textup{Im}\left\\{\sqrt{(ks)^{2}-(2\pi\ell-t)^{2}}\right\\}\geq 0$ when $t$
is real and for all $\ell$. In other words, the branches are selected such
that $\left.\sqrt{(ks)^{2}-(2\pi\ell-t)^{2}}\right|_{t=2\pi\ell}=ks$ for all
$\ell$.

We can further accelerate the convergence of (A.7) by using the following
asymptotic expansion as $\ell\rightarrow\infty$,

$\displaystyle\frac{2}{\sqrt{(ks)^{2}-(2\pi\ell-t)^{2}}}+\frac{2}{\sqrt{(ks)^{2}-(2\pi\ell+t)^{2}}}+\frac{2i}{\pi\ell}=\frac{2t^{2}+(ks)^{2}}{4i(\pi\ell)^{3}}+O(\ell^{-5}).$
(A.8)

Adding and subtracting the right hand side terms of (A.8) into (A.7) will give
us the accelerated formula,

$\displaystyle K(e^{it})=$ $\displaystyle\
H^{(1)}_{0}(ka)-1-\frac{2i}{\pi}\left(\gamma+\ln\left(\frac{ks}{4\pi}\right)\right)+\frac{2}{\sqrt{(ks)^{2}-t^{2}}}+\frac{2t^{2}+(ks)^{2}}{4i\pi^{3}}\zeta(3)$
$\displaystyle+\sum_{\ell=1}^{\infty}\Bigg{(}\frac{2}{\sqrt{(ks)^{2}-(2\pi\ell-t)^{2}}}+\frac{2}{\sqrt{(ks)^{2}-(2\pi\ell+t)^{2}}}+\frac{2i}{\pi\ell}-\frac{2t^{2}+(ks)^{2}}{4i(\pi\ell)^{3}}\Bigg{)},$
(A.9)

where $\zeta(z)$ is the Riemann zeta function. More terms in the expansion
(A.8) will allow us the accelerate the convergence even further,

$\displaystyle K(e^{it})=$ $\displaystyle\
H^{(1)}_{0}(ka)-1-\frac{2i}{\pi}\left(\gamma+\ln\left(\frac{ks}{4\pi}\right)\right)+\frac{2}{\sqrt{(ks)^{2}-t^{2}}}+\frac{2t^{2}+(ks)^{2}}{4i\pi^{3}}\zeta(3)+\frac{8t^{4}+24(kst)^{2}+3(ks)^{4}}{64i\pi^{5}}\zeta(5)$
$\displaystyle+\sum_{\ell=1}^{\infty}\Bigg{(}\frac{2}{\sqrt{(ks)^{2}-(2\pi\ell-t)^{2}}}+\frac{2}{\sqrt{(ks)^{2}-(2\pi\ell+t)^{2}}}+\frac{2i}{\pi\ell}-\frac{2t^{2}+(ks)^{2}}{4i(\pi\ell)^{3}}-\frac{8t^{4}+24(kst)^{2}+3(ks)^{4}}{64i(\pi\ell)^{5}}\Bigg{)}.$
(A.10)

If we take (A.7) with $N$ extra corrections and truncated at $L$, the
neglected tail-end of that sum has the asymptotic behaviour
$O\left(L^{-(2N+2)}\right)$ as $L\rightarrow\infty$, which is also the rate of
convergence. Figure 8 additionally compares (A.7), (A.9) and (A.10) against
the method of tail-end asymptotics. It shows that these formula are generally
much faster to converge than (A.5), until it gets to machine precision and
then the working precision needs to be increased.

Figure 8: These graphs display the progression of the different Wiener–Hopf
kernel methods (left) at $z=1$ and their absolute error with a high order
reference value (right) with respect to the truncation point $L$. Here we have
$ks=1$ and $ka=0.01$. For both figures $K(1)$ is computed using either the
original formula (A.1) (see blue line), with $N$ additional tail-end
corrections (A.5) ($N=1,2$ and $4$ for the red, green and orange lines
respectively) or using the fast-convergent formula ((A.7) is the black line,
(A.9) is magenta and (A.10) is cyan). The legend also gives the rate of
convergence with respect to $L$ for each line.

Another benefit of (A.7) is the ability to converge for complex values of $t$
(i.e. an annulus about the $z$ unit circle). This ability coupled with the
fast convergence, is why we use (A.10) for analytic evaluations of the
Wiener–Hopf kernel.

###  A.3 Rational approximation

To produce the final results of this article, we will need to create a
rational approximation of the discrete Wiener–Hopf kernel. After construction,
the rational approximation is computationally very fast, remarkably accurate
and is also trivial to factorise analytically for the Wiener–Hopf technique.
This rational approximation is given by

$\displaystyle
K(z)\approx\widetilde{K}(z)=\widetilde{K}_{1}\prod_{\ell=1}^{L}\frac{(z-z_{\ell})}{(z-p_{\ell})},$
(A.11)

where the constant $\widetilde{K}_{1}$ and the zero and pole locations
($z_{\ell}$ and $p_{\ell}$ respectively) are determined for the most accuracy.
In MATLAB, the function we use to create this rational approximation is “AAA”
in Chebfun. Here, we input the values of $K(z)$ on a set of points located on
the $z$ unit circle using (A.10) and then the “AAA” function outputs the
rational approximation of the form (A.11) as well as individual lists of zeros
($z_{\ell}$), poles ($p_{\ell}$) and associated residues.

Figure 9 compares the exact formula (A.10) and the “AAA” approximation (A.11)
on the $z$ complex plane by plotting the complex argument Wegert [2012]. This
comparison shows that the zero and poles of the rational approximation are
positioned to best imitate the branch cuts.

Figure 9: Phase portrait comparison between the exact formula for the
Wiener–Hopf kernel (A.10) (left) and the “AAA” approximation (A.11) (right).
Phase portraits indicate the complex argument of a function by using a HSV
colour scale. In this example, we have $ka=0.01$ and $ks=1$.

For convenience, the zeros and poles of the approximation are split into two
sets depending on their locations being inside or outside the unit circle.
Hence, we split and relabel the lists $z_{\ell}$ and $p_{\ell}$ with
superscripts $-$ and $+$ to indicate their position inside and outside
$|z|=1$,

$\displaystyle
K(z)\approx\widetilde{K}(z)=\widetilde{K}_{1}\prod_{\ell=1}^{L^{-}}\frac{(z-z^{-}_{\ell})}{(z-p^{-}_{\ell})}\prod_{\ell=1}^{L^{+}}\frac{(z-z^{+}_{\ell})}{(z-p^{+}_{\ell})},$
(A.12)

where $L^{-}+L^{+}=L$. In these lists, we expect
$\arg(z^{+}_{\ell})\approx\arg(p^{+}_{\ell})\approx\arg(e^{-iks})$ and
$\arg(z^{-}_{\ell})\approx\arg(p^{-}_{\ell})\approx\arg(e^{iks})$. We also
expect that $|z^{+}_{\ell}|,|p^{+}_{\ell}|>1$ and
$|z^{-}_{\ell}|,|p^{-}_{\ell}|<1$ and then, without loss of generality, the
zeros and poles are organised by absolute value from the unit circle, i.e.

$\displaystyle|z^{+}_{\ell}|<|z^{+}_{\ell+1}|,\ \
|p^{+}_{\ell}|<|p^{+}_{\ell+1}|,\ \ |z^{-}_{\ell}|>|z^{-}_{\ell+1}|,\ \
|p^{-}_{\ell}|>|p^{-}_{\ell+1}|.$ (A.13)

For this approximation to be effective, it is very important to try to
conserve the kernel identity $K(z)=K\left(\frac{1}{z}\right)$ as much as
possible. For this, we need the same number of zeros and poles inside the unit
circle as outside (meaning $L^{-}=L^{+}=L/2$). Ideally, the zeros and poles
inside and outside should also be in reciprocal pairs (or at least
approximately so), in other words $z^{+}_{\ell}z^{-}_{\ell}\approx 1$ and
$p^{+}_{\ell}p^{-}_{\ell}\approx 1$ for all $\ell$ (see Figure 10 (left) for a
test case). We are not aware of a way to enforce these properties in the “AAA”
approximation, but we can increase the likelihood by using reciprocal pairs in
the $z$ input data and employing checks that restart the approximation with
fresh data when the properties are not satisfied enough. Figure 10 (right)
plots the absolute error of the approximation. This figure illustrates the
high accuracy rating of the rational approximation not only on the unit circle
from which it is created, but in its neighbourhood too (the exception being in
the vicinity of the branch points and cuts).

Figure 10: On the left side, we plot the absolute value of the two quantities
$z^{+}_{\ell}z^{-}_{\ell}\approx 1$ and $p^{+}_{\ell}p^{-}_{\ell}\approx 1$
for all $\ell$. On the right side, we plot the absolute error of the rational
approximation $|K(z)-\widetilde{K}(z)|$ where the colour map is scaled
logarithmically. In this example, we have $ka=0.01$ and $ks=1$.

The downside to this approximation is that it cannot capture the correct
asymptotic behaviour at the branch points, no matter how good the
approximation is elsewhere. However, this is to be expected considering the
format of the rational approximation.

When factorising the Wiener–Hopf kernel, we recall how the $+$ (or $-$) factor
must be analytic and zero-free inside (or outside) the unit circle. For the
approximate kernel, it is clear that the $+$ factor will have zeros at
$z^{+}_{\ell}$ and poles at $p^{+}_{\ell}$. Similarly the $-$ factor will have
zeros and poles at $z^{-}_{\ell}$ and $p^{-}_{\ell}$ respectively. The result
of the factorisation is given by,

$\displaystyle
K^{+}(z)\approx\widetilde{K}^{+}(z)=\widetilde{K}^{+}_{1}\prod_{\ell=1}^{L/2}\frac{z-z^{+}_{\ell}}{z-p^{+}_{\ell}},\quad
K^{-}(z)\approx\widetilde{K}^{-}(z)=\widetilde{K}^{-}_{1}\prod_{\ell=1}^{L/2}\frac{z-z^{-}_{\ell}}{z-p^{-}_{\ell}},$
(A.14)

where the two constants $\widetilde{K}^{+}_{1}$ and $\widetilde{K}^{-}_{1}$
satisfy $\widetilde{K}^{+}_{1}\widetilde{K}^{-}_{1}=\widetilde{K}_{1}$ by
design. We chose $\widetilde{K}^{+}_{1}$ such that
$\widetilde{K}^{+}(z)=\widetilde{K}^{-}(1/z)$ is satisfied on at least one
point on the unit circle. For example, say $z=1$ (or $z=-1$ to avoid errors
due to branch point singularities if $ks\approx 2\pi M$ ($M\in{\mathbb{Z}}$)),
which implies that,

$\displaystyle\widetilde{K}^{+}_{1}=\left(\widetilde{K}_{1}\prod_{\ell=1}^{L/2}\frac{(1\mp
z^{-}_{\ell})(1\mp p^{+}_{\ell})}{(1\mp p^{-}_{\ell})(1\mp
z^{+}_{\ell})}\right)^{\frac{1}{2}}=\left(\widetilde{K}(\pm
1)\right)^{\frac{1}{2}}\prod_{\ell=1}^{L/2}\frac{1\mp p^{+}_{\ell}}{1\mp
z^{+}_{\ell}},$ (A.15)

then $\widetilde{K}^{-}_{1}=\widetilde{K}_{1}/\widetilde{K}^{+}_{1}$. Note
that the upper/lower signs in (A.15) are the result of choosing to evaluate
the factors at $z=1$ and $-1$ respectively. Ideally, the optimal choice for
these constants would require the quantity
$\widetilde{K}^{+}(z)-\widetilde{K}^{-}(1/z)$ to be minimised for all $z$ on
the unit circle. However, experience with test cases have shown that
optimisation makes a negligible difference to the value of the constants.

Lastly, we can use the rational approximation to accurately approximate the
$\lambda_{n}$ coefficients. We perform a partial fraction expansion on the
reciprocal of $\widetilde{K}^{+}(z)$,

$\displaystyle\frac{1}{\widetilde{K}^{+}(z)}$
$\displaystyle=\frac{1}{\widetilde{K}^{+}_{1}}\left[1+\sum_{\ell=1}^{L/2}\frac{\textrm{r}_{\ell}}{z-z^{+}_{\ell}}\right],\
\ \textrm{where}\ \
\textrm{r}_{m}=(z^{+}_{m}-p^{+}_{m})\prod_{\begin{subarray}{c}\ell=1\\\
\ell\neq
m\end{subarray}}^{L/2}\frac{z^{+}_{m}-p^{+}_{\ell}}{z^{+}_{m}-z^{+}_{\ell}},$
(A.16)

and then we substitute (A.16) into (2.23) to get these approximations for
$\lambda_{n}$

$\displaystyle\lambda_{0}$
$\displaystyle\approx\frac{1}{\widetilde{K}^{+}_{1}}\prod_{\ell=1}^{L/2}\frac{p^{+}_{\ell}}{z^{+}_{\ell}},\
\ \textrm{and}\ \
\lambda_{n}\approx-\frac{1}{\widetilde{K}^{+}_{1}}\sum_{\ell=1}^{L/2}\frac{\textrm{r}_{\ell}}{(z^{+}_{\ell})^{n+1}}.$
(A.17)

To test the accuracy of the $\lambda_{n}$ approximations, Figure 11 plots the
$\lambda_{n}$ coefficients from $n=0$ to $1000$ on the complex plane using
both the rational approximation (A.17) and the integral methods (2.26).

Figure 11: A comparison of the $\lambda_{n}$ coefficients $(0\leq n\leq 1000)$
on the complex plane by the approximation formulae (A.17) (blue) and the
integral formulae (2.26) (red). Note that the outlier is $\lambda_{0}$. Here
we have $ks=1$ and $ka=0.01$.

We find that the two graphs are on top of each other. This is clear that
evaluating the $\lambda_{n}$ coefficients via (A.17) is very accurate and we
also find that it is significantly faster computationally than using the
integral formula. However, it will only be as accurate as the rational
approximation itself and numerical error will eventually catch up as
$n\rightarrow\infty$ due to the denominator’s exponent.