# Trapped acoustic waves and raindrops: high-order accurate integral equation
method for localized excitation of a periodic staircase

Fruzsina J. Agocs<EMAIL_ADDRESS>Alex H. Barnett
<EMAIL_ADDRESS>
###### Abstract

We present a high-order boundary integral equation (BIE) method for the
frequency-domain acoustic scattering of a point source by a singly-periodic,
infinite, corrugated boundary. We apply it to the accurate numerical study of
acoustic radiation in the neighborhood of a sound-hard two-dimensional
staircase modeled after the El Castillo pyramid. Such staircases support
trapped waves which travel along the surface and decay exponentially away from
it. We use the array scanning method (Floquet–Bloch transform) to recover the
scattered field as an integral over the family of quasiperiodic solutions
parameterized by their on-surface wavenumber. Each such BIE solution requires
the quasiperiodic Green’s function, which we evaluate using an efficient
integral representation of lattice sum coefficients. We avoid the
singularities and branch cuts present in the array scanning integral by
complex contour deformation. For each frequency, this enables a solution
accurate to around 10 digits in a couple of seconds. We propose a residue
method to extract the limiting powers carried by trapped modes far from the
source. Finally, by computing the trapped mode dispersion relation, we use a
simple ray model to explain an observed acoustic “raindrop” effect (chirp-like
time-domain response).

††journal: Journal of Computational Physics

[1]organization = Center for Computational Mathematics, Flatiron Institute,
addressline = 162 Fifth Avenue, city = New York, state = NY, country = USA

##  1 Introduction

Periodic surface geometries have long been exploited to manipulate
electromagnetic and acoustic waves. Examples on small and large length-scales
include photonic crystals [1, 2], acoustic metamaterials [3, 4], diffraction
gratings [5], antennae [6, 7], anechoic chambers and materials [8], and
amphitheatres [9]. However, the accurate numerical solution of radiation near
such geometries faces several challenges. The domain—a perturbation of the
upper half-space—is unbounded vertically, thus its truncation must be
incorporate the correct upward-going radiation conditions. The possibility of
coupling to waves trapped along the surface, which in two dimensions (2D) do
not decay at all, means that ${\mathcal{O}}(1)$ reflection errors would result
by naive truncation to any finite number of unit cell periods. These trapped
waves may also cause resonances with high parameter sensitivity [10]. Having
non-periodic excitation breaks the periodicity of the problem, which makes
periodization impossible at first glance. However, it is possible to
deconstruct the single point-source solution into sets of quasiperiodic
solutions (the array scanning or Floquet–Bloch method). As in any wave
problem, high frequencies may demand a large discretization density [11].
Finally, the staircase model that we focus on here introduces corner
singularities that must be addressed in any high-order solver. Since solvers
are often part of a design optimization loop to tune material or shape
parameters, they must be robust and efficient. However, the above difficulties
mean that nonconvergent Rayleigh expansions [12], or ray asymptotic methods
[13] are often used in acoustic settings.

The purpose of this paper is twofold: 1) We present a high-order numerical
boundary integral (BIE) method for 2D acoustic scattering of a point source
from a singly-periodic geometry, in particular combining it with a high-order
accurate array scanning method via contour deformation. 2) We apply this to
the accurate numerical study of acoustic radiation in the neighborhood of a
sound-hard 2D staircase model, in order to understand time- and frequency-
domain phenomena due to nearby point-source excitation. For the time-domain we
compute the dispersion relation of trapped modes and use it to explain an
observed chirp effect [14]; for the frequency domain we propose a method to
extract the amplitude of the left- and right-going trapped modes, using
residues in the complexified on-surface wavenumber. Architecturally, solid
staircases with sound-hard surfaces are clearly very common; yet, we have not
found an accurate numerical study of acoustic trapping and guiding in their
vicinity. By combining high-order corner quadratures and lattice sums, we show
that a BIE can achieve around 10 digits of accuracy with solution times of a
couple of seconds per excitation frequency.

We consider a 2D model for linearized acoustics in a simply-connected region
$\Omega\subset\mathbb{R}^{2}$ containing a constant-density gas with constant
sound speed $c>0$ [15, 16], [17, Sec. 3.1]. This region lies above a connected
corrugated boundary $\partial\Omega$ extending with spatial period $d$ in the
$x_{1}$ direction, and is unbounded in the positive $x_{2}$ direction; see
fig. 1. The boundary has a maximum and minimum $x_{2}$ coordinate. Our
paradigm example will be the slope-$\pi/4$ right-angle staircase shown in the
figure, whose repeating element (unit cell) comprises two equal-length line
segments at an angle $\pi/2$ to each other. We write
$\mathbf{x}=(x_{1},x_{2})$. This models a 3D situation in which both geometry
and source are invariant in the 3rd (out of plane) direction.

The tools we present mostly concern frequency domain solutions, but we will
also use them to understand a chirp phenomenon for the following time domain
problem. Consider an impulse excitation by a point source at time $t=0$ and
location $\mathbf{x}_{0}\in\Omega$; this is a good model for a clap or
footstep in the acoustic application. The acoustic pressure $U(t,\mathbf{x})$
then obeys the wave equation with a source,

$\frac{\partial^{2}U(t,\bm{\mathbf{x}})}{\partial t^{2}}-\Delta
U(t,\bm{\mathbf{x}})=\delta(t)\delta(\mathbf{x}-\mathbf{x}_{0}),\qquad
t\in\mathbb{R},\quad\mathbf{x}\in\Omega,$ (1)

where $\Delta=\nabla\cdot\nabla$ is the spatial Laplacian, and by a rescaling
of time $t$ we set the sound speed to $c=1$. We assume quiescence before the
excitation: $U\equiv 0$ for $t<0$. 111Note that 1 could thus be rephrased as
the homogeneous wave equation in $t>0$ with initial conditions
$U(0,\cdot)\equiv 0$, $U_{t}(0,\cdot)=\delta_{\mathbf{x}_{0}}$. The sound-hard
(Neumann) boundary condition is

$U_{n}(t,\mathbf{x})=0,\qquad
t\in\mathbb{R},\quad\mathbf{x}\in\partial\Omega,$ (2)

where $U_{n}\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{n}\cdot\nabla U$ is
the normal derivative, $\mathbf{n}$ being the unit boundary normal vector
pointing into $\Omega$. This arises physically since the normal component of
the fluid velocity vanishes; it is a good approximation to an air-solid
interface [18, Ch. 1].

Taking the Fourier transform of 1 and 2 with respect to $t$ gives the main
focus of this paper, the inhomogeneous Helmholtz Neumann boundary value
problem (BVP),

$\displaystyle-(\Delta+\omega^{2})u$ $\displaystyle=$
$\displaystyle\delta_{\mathbf{x}_{0}}\qquad\mbox{ in }\Omega,$ (3)
$\displaystyle u_{n}$ $\displaystyle=$ $\displaystyle 0\qquad\mbox{ on
}\partial\Omega,$ (4)

and radiation conditions explained below. Here $\omega\in\mathbb{R}$ is the
frequency, and the wavelength is $2\pi/\omega$. From the $\omega$-dependence
of its solution $u$ one may understand features of the above time-domain
solution $U$.

(a)

(b)

Figure 1: Left: Photograph of _El Castillo_ at Chichen Itza, Mexico. Center:
2D infinite staircase geometry with coordinates used. Right: a single period
$\Gamma$ of the boundary, now shown with $x_{1}$ horizontal (the orientation
used throughout).

The remainder of this extended Introduction sets up definitions needed to
convert 3 and 4 into BVPs posed on a bounded unit cell, and overviews the rest
of the paper. The solution to 3 and 4 may be decomposed into a linear
combination of _quasiperiodic_ solutions, via what is known as the array
scanning method in the engineering literature [6, 19, 20], or inverse
Floquet–Bloch transform in mathematics [21, 22]. This family of quasiperiodic
solution is parameterized by $\kappa\in\mathbb{R}$, the horizontal (along-
surface) wavenumber, and each member of the family has the phased
translational symmetry
$u_{\kappa}(x_{1}+nd,x_{2})=\alpha^{n}u_{\kappa}(x_{1},x_{2})$, for all
$n\in\mathbb{Z}$ and $\mathbf{x}\in\Omega$, where the _Bloch phase_ is
$\alpha=e^{i\kappa d}$. The set of _plane waves_ of the form
$e^{i\mathbf{k}\cdot\mathbf{x}}$ obeying the homogeneous Helmholtz equation
and quasiperiodic symmetry is discrete, namely
$\mathbf{k}_{n}=(\kappa_{n},k_{n})$, where

$\kappa_{n}\mathrel{\mathop{\mathchar 58\relax}}=\kappa+2\pi n/d,\qquad
n\in\mathbb{Z}$ (5)

is the shifted lattice of horizontal wavenumbers with the same $\alpha$, and

$k_{n}\mathrel{\mathop{\mathchar 58\relax}}=+\sqrt{\omega^{2}-\kappa_{n}^{2}}$
(6)

is the vertical wavenumber. $k_{n}$ is either positive real (upwards-
propagating), zero (horizontally propagating), or positive imaginary (upwards-
decaying).

Then, for a given $\kappa$, the quasiperiodic solution $u=u_{\kappa}$ solves
the BVP

$\displaystyle-(\Delta+\omega^{2})u$ $\displaystyle=f_{\kappa}\quad$
$\displaystyle\text{in }\Omega,\quad$ (PDE) (7) $\displaystyle u_{n}$
$\displaystyle=0\quad$ $\displaystyle\text{on }\partial\Omega,\quad$ (boundary
condition) (8) $\displaystyle u(x_{1}+d,x_{2})$ $\displaystyle=\alpha
u(x_{1},x_{2})$ $\displaystyle(x_{1},x_{2})\in\Omega,\quad$ (quasiperiodicity)
(9) $\displaystyle u(x_{1},x_{2})$
$\displaystyle=\sum_{n\in\mathbb{Z}}c_{n}e^{i(\kappa_{n}x_{1}+k_{n}x_{2})},\quad$
$\displaystyle x_{2}>x_{2}^{0},\quad$ (radiation condition) (10)

where $x_{2}^{(0)}$ is any height lying above the support of the source and
above $\partial\Omega$, and $c_{n}$ are amplitudes of the outgoing plane
waves. The source $f_{\kappa}$ must naturally also be quasiperiodic: it is a
quasiperiodized version of the desired right-hand side
$\delta(\mathbf{x}-\mathbf{x}_{0})$, namely

$f_{\kappa}(\mathbf{x})\mathrel{\mathop{\mathchar
58\relax}}=\sum_{n=-\infty}^{\infty}e^{in\kappa
d}\delta(\mathbf{x}-\mathbf{x}_{0}-n\mathbf{d}),$ (11)

where $\mathbf{d}\mathrel{\mathop{\mathchar 58\relax}}=(d,0)$ is the lattice
vector. Fixing $\omega>0$, the above BVP has a unique solution for all
$\kappa\in\mathbb{R}$ apart from possibly a discrete set (which will
correspond to trapped modes, and discussed in detail in section 3) [23]. The
above also extends to $\kappa\in\mathbb{C}$; see [10, Sec. 4.1].

Each member $u_{\kappa}$ of the above quasiperiodic solution family is the
acoustic response of the periodic geometry to an infinite array of phased
point sources. The key observation (giving “array scanning” its name) is that
one may cancel out all but the central source by integration over $\kappa$ in
the first Brillouin zone:

$\delta(\mathbf{x}-\mathbf{x}_{0})=\frac{d}{2\pi}\int_{-\pi/d}^{\pi/d}f_{\kappa}(\mathbf{x})\mathrm{d}\kappa.$
(12)

This exploits the fact that $\int_{-\pi/d}^{\pi/d}e^{in\kappa
d}\mathrm{d}\kappa=2\pi/d$ for $n=0$, and zero for $n\in\mathbb{Z},n\neq 0$.
Since this is the right-hand side in 7, then, by linearity, performing the
same integral over the solution family $u_{\kappa}$,

$u(\mathbf{x})=\frac{d}{2\pi}\int_{-\pi/d}^{\pi/d}u_{\kappa}(\mathbf{x})\mathrm{d}\kappa$
(13)

recovers $u$, the solution to 7 and 18 with upward-propagating radiation
conditions.

######  Remark 1.

In the case where $u_{\kappa}$ is not unique for some discrete real $\kappa$
values, the integral 13 has poles on the real axis corresponding to evanescent
waves trapped on the corrugated surface. General conditions for existence of
such Neumann boundary condition trapped modes are not known [24, p.11–12],
although they are well known, and proven to exist for certain geometries [25,
26]. The mathematical formulation of causal outgoing boundary conditions in
the presence of trapped modes is subtle (e.g. , see [27]); we define this by
the choice of which side of these poles the integration path must take [22].
This is presented in section 4.

The integration interval $\kappa\in[-\pi/d,\pi/d]$ is known as the _first
Brillouin zone_ , and covers the family of solutions. This is because
$u_{\kappa}$ is the same for all $\kappa$ in the set 5 (since $\alpha$ is the
same, and so is the set of plane waves in the radiation condition, up to
relabeling). The integration may also be thought of as over $\alpha$ on the
unit circle [22]. We present a _high-order accurate_ quadrature method for the
above integral 13. This is crucial for efficiency, because each such
quadrature node demands a new BVP solution of 7, 8, 9 and 10. Here, care is
needed due to the possibility of poles mentioned in the above Remark, but also
because the integrand contains two square-root singularities (with associated
branch cuts) at so-called _(Rayleigh–)Wood anomalies_ [28, 29, 30]. These
anomalies are defined as pairs $(\omega,\kappa)$ where $k_{n}=0$ for some
$n\in\mathbb{Z}$ in 6 (see fig. 2), resulting in a horizontal plane wave.
Although the quasiperiodic BVP remains well-posed, it poses challenges for
Green’s function based numerical methods [31, 32, 33, 34, 35]. Following Zhang
[22] we use contour deformation to complex $\kappa$ to avoid such anomalies
and poles, and also optimize the contour shape to minimize the number of nodes
needed; see section 4.

Boundary integral equation (BIE) methods [17, 36] are especially suited for
the high-order accurate solution of each quasiperiodic BVP. Since they operate
by first converting the PDE to a boundary integral, then discretizing it to
yield a dense linear system, they only require the evaluation of a 1D integral
instead of the discretization of a (truncated) infinite 2D domain. This
reduction of dimensionality significantly reduces the number of unknowns, and
allows for an easy increase in order of accuracy via high-order quadrature
rules. Furthermore, the staircase geometry has two corners per period, one of
which induces fractional-power-law singularities in $u$; these are easily
handled with BIE. In contrast, a finite differencing (FD) or finite elements
(FEM) require meshing of the domain in a manner respecting the corner
singularity, and explicit handling of the radiation condition 10. Both FD [21]
and FEM [22] have been combined with array scanning before, but the only
numerical method with high-order spatial accuracy was applied to a locally
perturbed periodic boundary [37].

There exist alternative boundary-based approaches. These include meshfree
methods such as the method of fundamental solutions (MFS) [38, 39, 40], and
the plane waves method [41]. These share common roots with BIE methods in that
the solution is constructed from Helmholtz solutions in the domain, but
generally give ill-conditioned systems. For scattering by periodic surfaces, a
family of methods exists based on the Rayleigh hypothesis (Rayleigh methods)
[42, p. 17], [43]. These assume that an expansion of the scattered field like
10 is valid close to and on the surface, and approximate the solution as a
truncation of 10. Fast approximate solutions are most commonly derived using
the Helmholtz–Kirchhoff approximation [44], which assumes that each point on
the boundary scatters as if it was a plane at a slope matching that of the
boundary; this represents a short-wavelength limit, and its validity has been
studied thoroughly [12].

To use boundary integral methods, one splits the physical solution as
$u=u_{i}+u_{s}$, where the _incident wave_ $u_{i}$ solves 7 with free-space
radiation conditions and no boundary conditions, and is thus known
analytically, whilst the unknown _scattered wave_ $u_{s}$ solves the
homogeneous version of 7, with inhomogeneous boundary condition
$(u_{s})_{n}=-(u_{i})_{n}$ on $\partial\Omega$, and 9 and 10. A BIE is then
used to solve this BVP for $u_{s}$, as presented in section 2.

A BIE formulation on the infinite surface $\partial\Omega$ would not be
numerically feasible, and its truncated solution converges very slowly since
the contributions of distant sources decay only like $1/\sqrt{r}$. However,
the computation of the quasiperiodic solution may be reduced to a single unit
cell of the boundary. This periodization process involves replacing in the
integral kernels the free-space Green’s function $\Phi(\mathbf{x})$, defined
as the radiative solution to

$-(\Delta+\omega^{2})\Phi(\mathbf{x})=\delta(\mathbf{x}),$ (14)

with the quasiperiodic Green’s function ${\Phi_{\mathrm{p}}}(\mathbf{x})$
defined as the upwards- and downwards-radiating solution to

$-(\Delta+\omega^{2}){\Phi_{\mathrm{p}}}(\mathbf{x})=\sum_{n=-\infty}^{\infty}\alpha^{n}\delta(x_{1}-nd)\delta(x_{2}).$
(15)

The periodic Green’s function is therefore an infinite phased sum of single
point-source Green’s functions. The sum is slowly convergent and cannot be
used directly [45]. Yet a wide range of methods exist for the rapidly-
convergent approximation of ${\Phi_{\mathrm{p}}}$, including reformulation in
terms of quickly convergent lattice sums [46]. In section 2.1, we discuss one
such efficient computation of ${\Phi_{\mathrm{p}}}$ based on [47]. By the use
of contour deformation to complex $\kappa$, we will avoid Wood anomalies where
15 does not exist.

With the periodic Green’s function in hand, section 2.2 describes the boundary
layer representation of the solution and the discretization of the boundary,
including the treatment of the singularities present in the boundary integral
at the corners of the boundary. In general, these may be handled analytically,
e.g. via Gauss–Jacobi quadrature [48], or conformal mapping [49]; or handled
to high order by generalized Gauss quadrature [50], recursive compressed
inverse preconditioning [51], or rational function approximation [52]. In
section 2.3 that follows, we describe how the total field is reconstructed
outside the unit cell, and verify the high-order accuracy of the solution with
convergence tests which exploit flux conservation.

######  Remark 2.

A distributed spatial source (right-hand side in 7) may also be handled by a
slight generalization of our framework: one numerically computes their
incident wave $u_{i}$ by convolution of the source function with
${\Phi_{\mathrm{p}}}$. The BIE solutions for $u_{s}$ then proceed as before.

In section 3 we will show how numerically to locate eigenparameters
$(\omega,\kappa)$ for which trapped modes exist. Such modes are
eigenfunctions, i.e., nontrivial homogeneous solutions to the BVP 7, 8, 9 and
10. As 1 implied, at each $\omega$, such $\kappa$ values are vital to know
since they induce poles in the array scanning (inverse Floquet–Bloch)
integral. Yet, the mode _dispersion relation_ —the dependence of $\omega$ vs
$\kappa$ for trapped modes—also will provide a model to understand the
“raindrop effect” (chirp-like response) for the time-domain problem 1 and 2.
For this reason, we present a strategy to solve for the trapped $\omega$ at a
given $\kappa$. This involves rootfinding on the Fredholm determinant of the
$\omega$-dependent integral operator, as done in [53]. A trapped mode may be
reconstructed from the eigenfunction of the Fredholm 2nd-kind integral
equation. From the dispersion relation, we derive their group velocities; the
speed at which a given trapped mode propagates along the surface. We use an
approximate ray model (neglecting amplitudes) to predict the arrival times of
different frequencies at the bottom of a staircase modeled after the _El
Castillo_ pyramid, in order to understand the chirp-like sounds observed.

We also pose, and answer in Section 5, the following: how can one efficiently
use the array scanning method to report the power carried away by (left- or
right-going) trapped modes, as opposed to power radiated upwards away from the
surface? Characterizing this division of radiated power as a function of
frequency is crucial in related engineering applications such as point-source
radiators in nanophotonics and acoustic metamaterials. Finally, we discuss
avenues for future work in section 6.

##  2 BIE formulation, periodization, and discretization

Here we present the numerical method for solving the quasiperiodic BVP 7, 8, 9
and 10, at a given $\omega>0$ and $\kappa\in\mathbb{C}$. The physical wave
(potential) is written $u=u_{i}+u_{s}$. Fixing the source
$\mathbf{x}_{0}\in\mathbb{R}^{2}$, we take the incident wave to be the
quasiperiodic function

$u_{i}(\mathbf{x})={\Phi_{\mathrm{p}}}(\mathbf{x},\mathbf{x}_{0}),\quad\mathbf{x}\in\mathbb{R}^{2},$
(16)

where we use the notation
${\Phi_{\mathrm{p}}}(\mathbf{x},\mathbf{y})={\Phi_{\mathrm{p}}}(\mathbf{x}-\mathbf{y})$,
recalling 15. The PDE and boundary condition for $u$ 7 and 8 imply that
$u_{s}$ solves the BVP

$\displaystyle(\Delta+\omega^{2})u_{s}$ $\displaystyle=0\quad$
$\displaystyle\text{in }\Omega,$ (17) $\displaystyle(u_{s})_{n}$
$\displaystyle=-(u_{i})_{n}\quad$ $\displaystyle\text{on }\partial\Omega,$
(18)

with $u_{s}$ obeying the quasiperiodicity and radiation conditions 9 and 10.
Since the PDE is now homogeneous, a BIE solution becomes possible.

The unknown scattered field $u_{s}$ is represented by a quasiperiodic single-
layer potential

$u_{s}(\mathbf{x})={\mathcal{S}}\sigma=\int_{\Gamma}{\Phi_{\mathrm{p}}}(\mathbf{x},\mathbf{y})\sigma(\mathbf{y})\mathrm{d}s_{\mathbf{y}},\quad\mathbf{x}\in\mathbb{R}^{2},$
(19)

where the usual fundamental solution $\Phi(\mathbf{x},\mathbf{y})$ has been
replaced with its quasiperiodic counterpart, and
$\Gamma\mathrel{\mathop{\mathchar 58\relax}}=\partial\Omega\cap\\{|x_{1}|\leq
d/2\\}$ is a single unit cell (period) of the boundary. As $\mathbf{x}$
approaches the boundary from either side, let us define

$\displaystyle u_{s}^{\pm}(\mathbf{x})$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}=\lim_{h\to
0^{\pm}}u_{s}(\mathbf{x}+h\bm{\mathbf{n}}_{\mathbf{x}}),\quad\mathbf{x}\in\Gamma,$
(20) $\displaystyle(u_{s})_{n}^{\pm}(\mathbf{x})$
$\displaystyle\mathrel{\mathop{\mathchar 58\relax}}=\lim_{h\to
0^{\pm}}\bm{\mathbf{n}}_{\mathbf{x}}\cdot\nabla
u_{s}(\mathbf{x}+h\bm{\mathbf{n}}_{\mathbf{x}}),\quad\mathbf{x}\in\Gamma.$
(21)

These then satisfy the jump relations (see [36, Ch 6.3])

$\displaystyle u_{s}^{\pm}(\mathbf{x})$
$\displaystyle=S\sigma,\quad\mathbf{x}\in\Gamma,$ (22)
$\displaystyle(u_{s})_{n}^{\pm}(\mathbf{x})$
$\displaystyle=\left(D^{\mathrm{T}}\mp\frac{1}{2}I\right)\sigma,\quad\mathbf{x}\in\Gamma,$
(23)

where $D^{T}$ denotes the adjoint double-layer operator on $\Gamma$, namely
the operator with kernel
$\partial{\Phi_{\mathrm{p}}}(\mathbf{x},\mathbf{y})/\partial\bm{\mathbf{n}}_{\mathbf{x}}$
taken in the principal value sense. Here $\bm{\mathbf{n}}_{\mathbf{x}}$ is the
unit normal to the boundary at the target point $\mathbf{x}$.

The single-layer representation for $u_{s}$ automatically satisfies 17, and
the boundary condition can be written in terms of the single-layer jump
condition as

$(I-2D^{\mathrm{T}})\sigma=2(u_{i})_{n}\quad\text{on }\Gamma.$ (24)

This is a Fredholm integral equation of the second kind.

######  Remark 3.

The expert reader may wonder why a combined-field representation (“CFIE”) is
not needed here to prevent spurious resonances, as in the case of a bounded
obstacle [15]. In fact 19 is sufficient when the unbounded $\partial\Omega$ is
a graph of a function, because then the Dirichlet BVP in the complementary
domain $\mathbb{R}^{2}\backslash\overline{\Omega}$ is unique for any $\omega$;
see 4 below and its proof.

Figure 2: Location of Wood anomalies (black lines) and the first Brillouin
zone (red shaded region of the $\kappa$ axis) in the $\omega$-$\kappa$ plane,
for the case of spatial periodicity $d=1$.

###  2.1 Evaluation of quasiperiodic Green’s functions

Here we describe an efficient numerical scheme to evaluate
${\Phi_{\mathrm{p}}}(\mathbf{x})$ appearing in the above BIE, using a local
expansion about the origin, with lattice sum coefficients evaluated following
[47]. Note that, since it is only valid up to a vertical height comparable to
the unit cell width, a different plane-wave representation given in section
2.3 will be used for $x_{2}$ values beyond this. Recall the definitions 5 and
6, and that a Wood anomaly is a point in the $(\omega,\kappa)$ parameter plane
where $\kappa_{n}$ crosses $\kappa_{n}=\omega$. At such anomalies
${\Phi_{\mathrm{p}}}$ does not exist (e.g. see [45, 33]). We will thus assume
that $(\omega,\kappa)$ is not at a Wood anomaly (and direct the reader to [33,
34] for methods for quasiperiodic problems precisely at Wood anomalies).

The quasiperiodic Helmholtz Green’s function in 2D is given by [45]

${\Phi_{\mathrm{p}}}(x_{1},x_{2};\omega,d,\alpha_{0})=\frac{i}{4}\sum_{n=-\infty}^{\infty}H_{0}^{(1)}\left(\omega\sqrt{(x_{1}-nd)^{2}+x_{2}^{2}}\right)\alpha^{n}.$
(25)

The sum’s slow, conditional convergence means it is of little practical use.
We follow the approach described in [47] to rewrite it in terms of rapidly
convergent lattice sums. Our derivation, however, differs in two key ways; we
therefore reproduce some of the calculation below for convenience. Graf’s
addition theorem [54, (10.23.7)] is used first to expand each of the
$H_{0}^{(1)}$ centered on source point with an $x_{1}$ coordinate outside of
the unit cell, $|x_{1}|>\tfrac{d}{2}$, around the equivalent point inside the
unit cell. The resulting expression contains the Bessel functions
$H_{n}^{(1)}$ and $J_{n}$. Terms multiplying the same-order $J_{n}$ are
collected; their coefficient, the $n$th order lattice sum $S_{n}$, is
expressed as a contour integral. For faster convergence, following [33], we
instead exclude $|n|\leq 1$ from the sum in 25 and add them together directly;
this can be thought of splitting the periodic Green’s function into a near and
far component:

${\Phi_{\mathrm{p}}}={\Phi_{\mathrm{p},\text{near}}}+{\Phi_{\mathrm{p},\text{far}}},$
(26)

with

${\Phi_{\mathrm{p},\text{near}}}=\frac{i}{4}\left[H_{0}^{(1)}\left(\omega\rho_{0}\right)+\alpha^{-1}H_{0}^{(1)}\left(\omega\rho_{1}\right)\
+\alpha H_{0}^{(1)}\left(\omega\rho_{-1}\right)\right],$ (27)

where $\rho_{j}=\sqrt{(x_{1}-jd)^{2}+x_{2}^{2}}$. We use Graf’s theorem on the
remaining terms to write

${\Phi_{\mathrm{p},\text{far}}}=\frac{i}{4}S_{0}(\omega
d,\kappa)J_{0}(\omega\rho_{0})+\frac{i}{2}\sum_{n=1}^{\infty}S_{n}(\omega
d,\kappa)J_{n}(\omega\rho_{0})\cos(n\phi),$ (28)

where $\phi$ is the angle describing $(x_{1},x_{2})$ in polar coordinates. The
$S_{n}$ are the lattice sums of $n$th order, they only depend on the geometry
of the boundary and not the location of the target, therefore they only need
to be computed once for each value of $\omega$ and $\kappa$. We follow [47] to
derive the integral representation

$\displaystyle S_{n}(\omega d,\kappa)\approx-
ie^{i\tfrac{\pi}{4}}\frac{\sqrt{2}}{\pi}\bigg{[}$
$\displaystyle(-1)^{n}\alpha^{-2}\int_{0}^{a}[G_{n}(\tau)+G_{n}(-\tau)]F(\tau;\omega
d,\kappa)\mathrm{d}t$
$\displaystyle+\alpha^{2}\int_{0}^{a}[G_{n}(\tau)+G_{n}(-\tau)]F(\tau;\omega
d,-\kappa)\mathrm{d}t\bigg{]},$ (29)

with $\tau=(1-i)t$ and

$\displaystyle G_{n}(t)$ $\displaystyle=\left(t-i\sqrt{1-t^{2}}\right)^{n},$
(30) $\displaystyle F(t;\omega d,\kappa)$ $\displaystyle=\frac{e^{2i\omega
d\sqrt{1-t^{2}}}}{\sqrt{1-t^{2}}[1-e^{i\omega d(\sqrt{1-t^{2}}-\kappa/(\omega
d))}]}.$ (31)

The integrals appearing in $S_{n}$ may be evaluated numerically using the
trapezoidal rule. The reason behind 29 being an approximate expression instead
of an equality is that upper limit of the integrals, $a$, has been truncated
from $\infty$, exploiting the fact that the integrand decays quickly for $t\gg
1$ due to the exponential factor in $F$.

Three convergence parameters are needed here: the number of terms $N$ at which
the sum in ${\Phi_{\mathrm{p},\text{far}}}$ (in 28) is truncated, the upper
limit $a$ of the integrals in 29, and the number of nodes $M$ used in the
trapezoidal rule to calculate 29. For the typical $\omega$ and accuracies we
present, we found $N=40$, $a=15$, and $M=10^{4}$ sufficient.

###  2.2 Discretizing the boundary integral

We discretize and solve the integral equation 24 using the Nyström method [36,
Sec 12.2], summarized below. First, 24 is written in the standard form

$\sigma(t)-\int_{a}^{b}K(t,s)\sigma(s)\mathrm{d}s=f(t)\quad\text{for}\quad
t\in[a,b],$ (32)

where $K(t,s)$ is the appropriate kernel function, and the boundary is
parameterized by $s$ that runs from $s=a$ to $s=b$. Then, the unknown
_density_ function $\sigma(t)$ is approximated by another function
$\sigma_{n}(t)$ that obeys

$\sigma_{n}(t)-\sum_{j=1}^{n}w_{j}K(t,s_{j})\sigma_{n}(s_{j})=f(t),$ (33)

where the integral has been replaced with a quadrature formula with nodes
$\\{s_{i}\\}_{i=1}^{n}$ and weights $\\{w_{i}\\}_{i=1}^{n}$, and the kernel
function $K$ with a rank-$n$ operator. Then the values of $\sigma_{n}(t)$ at
the nodes $\\{s_{i}\\}$, $\\{\sigma_{i}^{(n)}\\}\mathrel{\mathop{\mathchar
58\relax}}=\\{\sigma_{n}(s_{i})\\}_{i=1}^{n}$ satisfy the linear system

$\sigma_{i}^{(n)}-\sum_{j=1}^{n}w_{j}K(s_{i},s_{j})\sigma_{j}^{(n)}=f(s_{i}),\quad\forall
i=1,2,\ldots,n,$ (34)

or in matrix notation,

$A\bm{\mathbf{\sigma}}^{(n)}=\bm{\mathbf{f}},$ (35)

with $A_{ij}=\delta_{ij}-K(s_{i},s_{j})w_{j}$. Then if any vector
$\bm{\mathbf{\sigma}}^{(n)}=\\{\sigma_{i}^{(n)}\\}$ solves the above system,
then $\sigma^{(n)}(t)$ at any $t$ can be reconstructed as

$\sigma^{(n)}(t)=f(t)+\sum_{j=1}^{n}w_{j}K(t,s_{i})\sigma_{j}^{(n)}.$ (36)

Substituting the kernel associated with our exterior Helmholtz problem with
boundary $\partial\Omega\mathrel{\mathop{\mathchar 58\relax}}=z(s)$, the $A$
matrix in 35 becomes:

$A_{ij}=\delta_{ij}-2w_{j}z^{\prime}(s_{j})\frac{\partial{\Phi_{\mathrm{p}}}(s_{i},s_{j})}{\partial\bm{\mathbf{n}}_{s_{i}}}\quad\forall
i,j=1,2,\ldots,n,$ (37)

while the components of $\bm{\mathbf{f}}$ are given by

$f_{j}=2(u_{i})_{n}(s_{j})\quad\forall j=1,2,\ldots,n.$ (38)

Finally, the gradient of ${\Phi_{\mathrm{p}}}$, necessary for filling the
$A$-matrix in 37, is

$\displaystyle\nabla{\Phi_{\mathrm{p}}}(\bm{\mathbf{x}},\bm{\mathbf{y}})=$
$\displaystyle\frac{i}{4}\bigg{(}-\omega
H_{1}^{(1)}(\omega\rho_{0})\frac{\bm{\mathbf{x}}-\bm{\mathbf{y}}}{\rho_{0}}-\omega
H_{1}^{(1)}(\omega\rho_{1})\frac{\bm{\mathbf{x}}-\bm{\mathbf{y}}-d\bm{\mathbf{e}}_{x_{1}}}{\rho_{1}}$
$\displaystyle-\omega
H_{1}^{(1)}(\omega\rho_{-1})\frac{\bm{\mathbf{x}}-\bm{\mathbf{y}}+d\bm{\mathbf{e}}_{x_{1}}}{\rho_{-1}}-\omega
J_{1}(\omega\rho_{0})S_{0}\frac{\bm{\mathbf{x}}-\bm{\mathbf{y}}}{\rho_{0}}$
(39) $\displaystyle+2\sum_{n=1}^{\infty}\omega
J_{n}^{\prime}(\omega\rho_{0})\cos(n\phi)S_{n}\frac{\bm{\mathbf{x}}-\bm{\mathbf{y}}}{\rho_{0}}+nJ_{n}(\omega\rho_{0})\sin(n\phi)S_{n}\frac{\hat{\bm{\mathbf{n}}}_{\bm{\mathbf{x}}-\bm{\mathbf{y}}}}{\rho_{0}}\bigg{)},$

where $\rho_{j}\mathrel{\mathop{\mathchar
58\relax}}=\sqrt{(x_{1}-y_{1}-jd)^{2}+(x_{2}-y_{2})^{2}}$,
$\hat{\bm{\mathbf{n}}}_{\bm{\mathbf{x}}-\bm{\mathbf{y}}}\mathrel{\mathop{\mathchar
58\relax}}=(x_{2}-y_{2},-(x_{1}-y_{1}))/\rho_{0}$ is a unit normal to the
difference vector $\bm{\mathbf{x}}-\bm{\mathbf{y}}$, and
$\bm{\mathbf{e}}_{x_{1}}$ is a unit vector in the $x_{1}$ direction.

The key task for achieving high-order accuracy is in choosing the appropriate
quadrature nodes and weights appearing in 33. For smooth boundaries, efficient
quadrature rules are known, based on either the global periodic trapezoid
rule, or high-order panel-wise quadrature. For boundaries with corners,
however, the density function $\sigma$ is singular at the corners. The
closeness of this singularity limits the radius of the Bernstein ellipse
associated with the interpolant used in Gaussian quadrature rules, which in
turn degrades the accuracy of interpolation and quadrature [55, Ch 8, 19]. To
maintain the order of accuracy whilst keeping the number of quadrature nodes
the same, the size of the quadrature interval needs to be reduced as the
corned is approached, i.e. the quadrature grid needs to be refined according
to the features of the boundary [56]. We divide the two sides of unit cell
boundary into equally-spaced panels, and use $P$-th order Gauss–Legendre nodes
and weights on each, i.e. $s_{i}$ and $w_{i}$ are always the Gauss–Legendre
nodes and weights on the standard interval $[-1,1]$. Then for each panel with
a corner as an endpoint, we divide the panel in a $1\mathrel{\mathop{\mathchar
58\relax}}(r-1)$ ratio with $r\geq 2$, using the same $P$-th order
Gauss–Legendre scheme on each. With $r=2$ (dyadic refinement), the net error
after $\mathcal{O}(\log_{2}(1/\varepsilon))$ levels of refinement will be
$\mathcal{O}(e^{-P}\log_{2}(1/\varepsilon)$, where $\varepsilon$ is a given
numerical precision [56]. Moving from a single set of $n$ global quadrature
nodes to panel quadrature with $Q$ panels, each with $P$ nodes, modifies the
expression 34 to

$\sigma_{j,i}-\sum_{q=1}^{Q}\sum_{p=1}^{P}w_{j,i,q,p}K(s_{j,i},s_{q,p})\sigma_{q,p}=f_{j,i},\quad\forall
j=1,\ldots,Q,\;i=1,\ldots,P,$ (40)

where $s_{j,i}$ is the $i$-th Gauss–Legendre node on the $j$-th panel,
$\sigma_{j,i}$, $f_{j,i}$, $w_{j,i}$ are the corresponding density, boundary
data, and quadrature weight, and $K(s_{j,i},s_{q,p})$ is the value of the
kernel associated with the given pair of quadrature nodes.

Typically, different quadrature rules are needed depending on where $s_{j,i}$
and $s_{q,p}$ (the “target” and “source” points) lie relative to each other
[50, 57], with special care required if they lie on adjacent or the same
panel. The need for this can be seen by inspecting 39 as
$\bm{\mathbf{x}}\to\bm{\mathbf{y}}$: simple Gauss–Legendre rules cannot
capture the (weak) singularity that emerges in this limit, therefore the
close-to-diagonal entries in $A$ will not be accurate. In the special case of
the staircase, however, it is possible to achieve higher order accuracy
without invoking special quadrature rules: since quadrature panels are linear,
the interaction between source and target points on the same panel or on the
same side of the boundary, is zero. This is again seen from 37 and 39:
$\bm{\mathbf{x}}-\bm{\mathbf{y}}$ is perpendicular to
$\bm{\mathbf{n}}_{\mathbf{x}}$ if $\bm{\mathbf{x}}$ and $\bm{\mathbf{y}}$ are
colinear.

Apart from the closeness of singularities, decrease in accuracy can arise from
catastrophic cancellation if
$|\bm{\mathbf{x}}-\bm{\mathbf{y}}|\ll|\bm{\mathbf{x}}|,|\bm{\mathbf{y}}|$,
e.g. if they point to quadrature nodes on small panels on either side of a
corner. To mitigate this, we parameterize the quadrature panels such that the
quadrature nodes’ positions are measured relative to the nearest corner.

Note that reconstructing the solution according to 36 close to the boundary
(i.e. with $t\to s_{i}$) would require special quadrature rules for close
evaluation since the position vectors of the target and source points are not
necessarily colinear. Since we do not employ any special rules, the numerical
solution will not be reliable within around the width of the smallest
quadrature panel of the boundary.

Figure 3: Discretization nodes on a single unit cell $\Gamma$ of the boundary.
The underlying coarse discretization has $8$ equal panels on each straight
line. Panels touching corners have then been subdivided dyadically $10$ times,
with shrinkage ratio $r=2$. Each resulting panel was populated with $16$
Gauss–Legendre quadrature nodes. An inset shows the result of the refinement.

###  2.3 Reconstructing the solution

Once the density evaluated at the boundary nodes,
$\bm{\mathbf{\sigma}}^{(n)}$, is obtained, on may compute the scattered field
$u_{s}$ at a point $\mathbf{x}$ inside the unit cell from its single-layer
representation,

$u_{s}(\mathbf{x})=\sum_{j=1}^{n}w_{j}z^{\prime}(s_{j}){\Phi_{\mathrm{p}}}(\mathbf{x},s_{j})\sigma^{(n)}_{j}.$
(41)

(After which, $u(t)=u_{i}(t)+u_{s}(t)$.) The lattice sum representation of
${\Phi_{\mathrm{p}}}$, and hence the above expression, quickly loses accuracy
outside of the circle $(x_{1}+x_{2})^{2}=\tfrac{3d}{2}$ due to the application
of Graf’s addition theorem222This would be a smaller, radius-$\tfrac{d}{2}$
circle were $n=\pm 1$ not excluded from 28. .

Above the unit cell, the upwards propagating radiation condition may be
exploited as follows. Let $x_{j}\mathrel{\mathop{\mathchar
58\relax}}=-\tfrac{d}{2}+\tfrac{jd}{N}$ for $j=0,1,\ldots,N-1$ and
$u^{0}_{j}\mathrel{\mathop{\mathchar 58\relax}}=u(x_{j},x_{2}^{0})$.
Evaluating the radiation condition at $x_{2}=x_{2}^{0}$, notice that after
rearranging, it takes the form of a discrete Fourier transform:

$\displaystyle u^{0}_{j}e^{-i\kappa
x_{j}}=\sum_{n\in\mathbb{Z}}\tilde{c}_{n}e^{i2\pi nx_{j}},$ (42)

with $\tilde{c}_{n}=c_{n}e^{ik_{n}x_{2}^{0}}$. Therefore, using the convention
from [58, Ch. 12.1],

$\tilde{c}_{n}=\frac{1}{N}\mathcal{F}\left(\\{u^{0}_{j}e^{-i\kappa
x_{j}}\\}_{i=0}^{N-1}\right)_{n},$ (43)

and

$u(x_{j},x_{2})=N\cdot\mathcal{F}^{-1}\left(\\{\tilde{c}_{n}e^{i(k_{n}(x_{2}-x_{2}^{0})+i\kappa
x_{j})}\\}_{n=0}^{N-1}\right)_{j}.$ (44)

In the $n$th neighboring unit cell, i.e.
$\tfrac{(2n-1)d}{2}<x_{1}\leq\tfrac{(2n+1)d}{2}$ (at any given $x_{2}$), the
solution is found using the quasiperiodicity condition,

$u(x_{1}+nd,x_{2})=\alpha^{n}u(x_{1},x_{2}).$ (45)

###  2.4 Convergence check via flux conservation

In the absence of sources of sinks inside a closed boundary $\Gamma$, the net
flux leaving $\Gamma$ is zero. How close it is to zero numerically can be used
to measure the accuracy of the method and investigate convergence in terms of
the size of the linear system 35. The net flux through the surface is

$F_{\Gamma}=\imaginary\left(\int_{\Gamma}\bar{u}u_{n}\mathrm{d}s\right)=0.$
(46)

To test the convergence of the method above, we choose the incident wave to be
a plane wave, $u_{i}=e^{i\omega\bm{\mathbf{n}}_{i}\cdot\bm{\mathbf{x}}}$,
where $\bm{\mathbf{n}}_{i}$ is a unit vector describing its direction of
travel. This corresponds to the limit of moving a quasiperiodic point source
array to infinity in the vertical direction. Specifically,
$\bm{\mathbf{n}}_{i}=(\cos(\phi_{i}),\sin(\phi_{i}))$ with
$-\pi\leq\phi_{i}\leq 0$. Then $\kappa=\omega\cos(\phi_{i})$. Let $\Gamma$ be
the perimeter of the unit cell as shown in fig. 1, bounded from below by the
stair boundary, above by the line $x_{2}=\tfrac{d}{2}$, and from the two sides
by $x_{1}=\pm\tfrac{d}{2}$. The net flux through the two sides has to be zero
by symmetry, and it is also zero across the lower edge of $\Gamma$ due to
Neumann boundary conditions. Therefore we have

$F_{\Gamma}=\imaginary\left(\int_{-d/2}^{d/2}\bar{u}\frac{\partial u}{\partial
x_{2}}\mathrm{d}x_{1}\right)=0.$ (47)

In fig. 4 we evaluate $F_{\Gamma}$ as we refine the corner-adjacent panels
with $r=3$. Initially there are $Q=8$ panels on each side with $P=8$ nodes on
each, and after $R$ levels of refinement, there are $Q=2(8+2R)$ panels and a
total of $PQ=128+32R$ nodes. The experiment used $\omega=4.0$, and
$\phi_{i}=-0.5$. The figure confirms exponential convergence and suggests that
around $10$ levels of refinement are needed for $\approx 7$ digits of
accuracy, and around $22$ for $10$ digits, at which point the system size is
still smaller than $10^{3}\times 10^{3}$. The reduction of error stops at this
point, which is roughly consistent with the spacing between the closest
quadrature nodes approaching machine precision.

Figure 4: Convergence test via flux conservation for the case of an incident
plane wave. The figure shows the net flux exiting the central unit cell, which
is analytically zero, as a function of both the number of times the corner-
adjacent quadrature panel has been subdivided on the boundary, and the size of
the resulting linear system. During refinement the panels have been split in a
$1\mathrel{\mathop{\mathchar 58\relax}}2$ ratio ($r=3$).

(a)  (b)

Figure 5: Left: numerically computed band structure (dispersion relation) for
evanescent trapped modes of the $\pi/4$-slope staircase with period $d=1$.
Only the right (positive) half of the Brillouin zone is shown. Dots show the
band structure $\omega_{\mathrm{tr}}(\kappa)$, and the line shows the values
$\omega=\kappa$. Radiation into the upper half plane is only possible when
$\omega>|\kappa|$. Right: group velocity $v_{g}\mathrel{\mathop{\mathchar
58\relax}}=d\omega_{\mathrm{tr}}(\kappa)/d\kappa$ plotted over the same
domain.  Figure 6:  The real part of two example trapped modes $\phi$, for the
$\pi/4$ staircase with $d=1$. Left: the highest frequency mode, at
$\kappa=\pi$. Right: An intermediate frequency mode, at $\kappa=1.54$.  Figure
7: Graph of frequency (in Hz) vs time (in seconds since emission) observed a
distance $D$ along a sound-hard staircase from an impulsive source for the
wave equation. Parameters for the _El Castillo_ staircase are used, with $D$
the full 91 steps comprising one side of the pyramid.

##  3 Trapped acoustic modes and the raindrop effect

In this section we describe a Fredholm determinant method to find evanescent
modes trapped by the corrugated sound-hard interface. Such modes are
eigenfunctions, i.e., nontrivial solutions to the homogeneous quasiperiodic
BVP 7, 8, 9 and 10, for some $\kappa\in\mathbb{R}$ and $\omega>0$. Their
parameters form a continuous families, thus curves in the $(\omega,\kappa)$
plane, known as the _band structure_ [2, 10]. Each curve may be described by
its trapped mode frequency $\omega_{\mathrm{tr}}(\kappa)$, commonly referred
to as a dispersion relation. For staircase geometries we find empirically that
there is only a single such trapped frequency at each wavenumber
$0<|\kappa|\leq d/2$ in the Brillouin zone, as shown in fig. 5.

Recall that the quasiperiodic BVP has the BIE formulation 24. With homogeneous
boundary data this becomes the BIE

$(I-2D^{T})\sigma=0,$ (48)

with $D^{T}$ the adjoint double-layer operator from the previous section. One
might then hope that the condition for a trapped mode to exist is equivalent
to the existence of a nontrivial density $\sigma$ solving 48. We now show that
this is indeed so, thus that the proposed numerical method is robust (free of
spurious resonances).

######  Theorem 4.

Fix $\omega>0$ and $\kappa\in\mathbb{R}$. There is a trapped mode (i.e., a
nontrivial $\phi$ solving the homogeneous quasiperiodic BVP 7, 8, 9 and 10) if
and only if $\operatorname{dim}\operatorname{Nul}(I-2D^{T})>0$.

###### Proof.

Let $\sigma$ be a nontrivial solution to $(D^{T}-I/2)\sigma=0$, and let
$\phi={\mathcal{S}}\sigma$ throughout $\mathbb{R}^{2}$, recalling that
${\mathcal{S}}$ is the quasiperiodic single-layer potential on $\Gamma$ (the
part of $\partial\Omega$ in the central unit cell). Then
$(\Delta+\omega^{2})\phi=0$ in $\mathbb{R}^{2}\backslash\partial\Omega$. By
the jump relations, $\phi_{n}^{+}=0$ (using the notations 20 and 21), and by
construction $\phi$ also satisfies 9 and 10. It remains to show that $\phi$ is
nontrivial. Assume this were so, then $\phi^{+}=0$, so by the jump relations,
$\phi^{-}=0$. However, $\phi$ in the half-space below $\partial\Omega$ would
then be a homogeneous solution to the (downward-facing) quasiperiodic
Dirichlet problem, which is unique [42, p. 56] [59, Thm. 2.1]. Thus $\phi$
would vanish below $\partial\Omega$, so $\phi_{n}^{-}=0$. By the jump relation
$\sigma=\phi_{n}^{-}-\phi_{n}^{+}$ would vanish, a contradiction with the
hypothesis. Thus $\phi$ is nontrivial, hence a trapped mode.

For the converse, let $\phi$ be a trapped mode, then by the quasiperiodic
version of Green’s representation theorem, in the upper domain
$\phi={\mathcal{S}}\phi_{n}^{+}+{\mathcal{D}}\phi^{+}$, where ${\mathcal{D}}$
is the quasiperiodic double-layer potential defined by

$({\mathcal{D}}\sigma)(\mathbf{x})\mathrel{\mathop{\mathchar
58\relax}}=\int_{\Gamma}\frac{\partial{\Phi_{\mathrm{p}}}(\mathbf{x},\mathbf{y})}{\partial\bm{\mathbf{n}}_{\mathbf{y}}}\sigma(\mathbf{y})\mathrm{d}s_{\mathbf{y}}~{}.$

Since $\phi_{n}^{+}=0$, and taking $\mathbf{x}$ to $\partial\Omega$ from
above, $\phi^{+}=(D+I/2)\phi^{+}$, showing that $D-I/2$ has a nontrivial null
vector. By the Fredholm–Riesz theory, the same holds for $D^{T}-I/2$, since it
is the adjoint with respect to the bilinear form
$\langle\psi,\phi\rangle\mathrel{\mathop{\mathchar
58\relax}}=\int_{\Gamma}\psi\phi\mathrm{d}s$ [17]. ∎

This informs our numerical approach: we fix $\kappa$ and solve a nonlinear
eigenvalue problem with respect to $\omega$. We adapt the Fredholm determinant
method of Zhao and the 2nd author [60]. $g(\omega)\mathrel{\mathop{\mathchar
58\relax}}=\det(I-2D^{T})$ is approximated by an $N\times N$ determinant,
where $I-2D^{T}$ is replaced by its Nyström matrix with entries given by 37.
At each $\kappa$, then $\omega_{\mathrm{tr}}(\kappa)$ is found as a root of
$g(\omega)$, using a simple Newton iteration with a convergence criterion
close to machine accuracy. For staircase geometries, we find only one such
root, and always with $\omega_{\mathrm{tr}}(\kappa)<|\kappa|$, implying that
the mode is trapped (frequency is below the light line), rather than embedded
in the continuous spectrum [2, 10]. Because the trapped frequencies never
intersect the light lines, the issue of Wood anomalies does not cause a
problem in the mode-finding task.

Figure 5 shows this set of $\omega_{\mathrm{tr}}$ found at each $\kappa$. The
gap between $\omega_{\mathrm{tr}}(\kappa)$ and $|\kappa|$ indicates the
strength of trapping (rapidity of evanescent decay as $x_{2}\to\infty$). The
most trapped mode is at $\kappa=\pm\pi$, known as an “optical mode”, and has
$\omega_{\mathrm{tr}}\approx 2.551$ for the $d=1$ staircase. As $\kappa\to 0$
we see trapping becoming arbitrarily weak.

With a zero of the determinant found, the null vector $\sigma$ obeying
$(I-2D^{T})\sigma=0$ is found via an SVD, then $\phi={\mathcal{S}}\sigma$
gives the mode. Numerically the lattice-sum expansion may only be used up to a
height around $d/2$ above the origin, so as before one must switch to Fourier
series evaluation above this. Figure 6 shows the real part of example modes.

###  3.1 Ray model for time-domain chirp response at El Castillo

The above trapped mode dispersion curve provides a simple explanation of the
acoustic “raindrop” effect on the long stone staircases at the _El Castillo_
pyramid at Chichen-Itza. Impulsive sources of sound, such as footsteps, are
reported to sound like short chirps (frequency rising vs time) when heard from
distances far up or down the staircase. Pending a full numerical investigation
of the time-domain solution to 1 and 2, we use a simple ray model for
dispersive propagation [61, Sec. 2.6]. Both source and receiver are assumed to
be close enough to the surface to couple well to the modes $\phi$ at all
wavenumbers (this is true for a footstep, less obviously so for a standing
listener). Absorption and modeling of amplitudes are ignored. Impulsive
excitation in 1 is assumed to excite all frequencies. Each frequency $\omega$
below the maximum of $\omega_{\mathrm{tr}}$ is partly trapped (see section 5),
and this component propagates along the staircase at its _group velocity_

$v_{g}=\frac{\mathrm{d}\omega_{\mathrm{tr}}}{\mathrm{d}\kappa},$ (49)

at the appropriate $\kappa_{\omega}$ such that
$\omega_{\mathrm{tr}}(\kappa_{\omega})=\omega$. We plot $v_{g}$ vs $\kappa$ on
the right of fig. 5, by differentiating the interpolant of
$\omega_{\mathrm{tr}}$ on the $\kappa$ grid. The frequency’s arrival time a
distance $D$ along the staircase is thus $t(\omega)=D/v_{g}(\omega)$. Finally,
inverting this last relationship gives the frequency $\omega$ heard at each
time $t$ after the impulse.

We now insert physical parameters for _El Castillo_. The nondimensionalized
speeds used in the rest of the paper must be multiplied the sound speed
$c=343$ m/s. According to [62] the stairs at this pyramid have depth equal to
height equal at $q=0.263$ m, i.e. a $\pi/4$ staircase with period
$d=\sqrt{2}q\approx 0.372$ m. The maximum trapped mode frequency is thus
$(c/d)2.551\approx 374$ Hz; this is the highest frequency explainable in the
model. 333The frequency of a free-space plane wave traveling along the
staircase with wavenumber at the Brillouin zone edge $\kappa=\pi/d$ is
$c/2d\approx 461$ Hz. Neither of these frequencies appear in [14], although a
“raindrop frequency” of 307.8 Hz is mentioned. Each staircase has 91 steps,
giving $D\approx 34$ m.

The resulting predicted frequency vs arrival time is shown in fig. 7. The
first arrivals are the lowest frequencies; for these the dispersion curve is
almost that of free air ($v_{g}\equiv c$) so they start arriving immediately
after any direct (non-trapped) radiation. Most of the frequency “chirp” occurs
during the first 0.2 s after first arrival, an interval containing of order 50
cycles, and thus plenty to detect the upwards frequency trend. A long “bell-
like” tail, asymptoting to the maximum 374 Hz, is expected, associated with
the slowly-propagating optical modes as shown on the left of fig. 6.

######  Remark 5.

We provide a link to a WAV file simulated to match the above chirp prediction
at this URL: https://doi.org/10.5281/zenodo.10005461. The authors do not have
access to audio recordings of footsteps on this staircase, and would welcome
the chance to validate the predictions.

##  4 Scattering from a single point source

As laid out in section 1, we simulate scattering from a single point source by
integrating over the solutions from a quasiperiodic array of point sources
with different $\kappa$ values. Let
$u_{\kappa,i}(\mathbf{x})={\Phi_{\mathrm{p}}}_{,\kappa}(\mathbf{x},\mathbf{x}_{0})$,
where the subscript $\kappa$ refers to having fixed $\kappa$ and $x_{0}$ is a
given source position within the central unit cell. Then from 11 and 12 it
follows that

$u_{i}(\mathbf{x})=\Phi(\mathbf{x},\mathbf{x}_{0})=\frac{1}{2\pi}\int_{-\pi}^{\pi}u_{\kappa,i}(\mathbf{x})\mathrm{d}\kappa,$
(50)

and with $u_{\kappa}(\mathbf{x})$ referring to the total field from an
incident field $u_{\kappa,i}(\mathbf{x})$ as above,

$u(\mathbf{x})=\frac{1}{2\pi}\int_{-\pi}^{\pi}u_{\kappa}(\mathbf{x})\mathrm{d}\kappa=\frac{1}{2\pi}\int_{-\pi}^{\pi}(u_{\kappa,i}(\mathbf{x})+u_{\kappa,s}(\mathbf{x}))\mathrm{d}\kappa,$
(51)

where $u_{\kappa,s}$ is the scattered field associated with $u_{\kappa,i}$.

If we were to evaluate the total field $u(\mathbf{x})$ at a point
$\mathbf{x}=(0.22,-0.16)$, fig. 8 shows the integrand that would appear in the
51 array scanning integral, in the complex $\kappa$ plane. The frequency was
set to $\omega=2.4$, and the single point source was located at
$x_{0}=(-0.2,0.1)$. Each side of the stair boundary was split into $4$ panels
initially, then refined $5$ times with $r=3$, until the total number of
quadrature panels over the boundary was $Q=28$, each with $P=16$ nodes. This
would give a result accurate to around $7$ digits according to the convergence
test in fig. 4. The two key features in the figure are the Wood anomalies
located at $\kappa=\pm\omega$ as branch points and two poles at
$\kappa=\pm\kappa_{\mathrm{tr}}$, with
$\omega\leq\kappa_{\mathrm{tr}}\leq\pi$. The branch cuts associated with the
branch points are indicated by dashed lines. The direction of the branch cuts
can be chosen by altering the integration path in 29 (by choosing $\tau(t)$),
but as it will become clear later, it is convenient to choose the cuts to lie
in the lower and upper half-plane for negative and positive $\kappa$,
respectively. If $\omega$ were larger than $\omega_{c}$, where no trapped
modes exist for any $\kappa$, there would be no poles present; and in the
limit of $\omega\to 0$ and $\omega\to\omega_{c}$, the two poles would coalesce
at $\kappa_{\mathrm{tr}}\to 0$ and $\kappa_{\mathrm{tr}}\to\pm\pi$,
respectively.

Figure 8: Real part of the integrand of the array scanning integral 51 plotted
in the complex $\kappa$-plane, where the field is evaluated at a single
target, $\mathbf{x}=(0.22,-0.16)$. Plotted over this are the branch points and
cuts (dots with dotted lines), poles (dots), and the contour deformation used
(solid black line).

The branch points and poles in 51 require careful treatment. If
$\omega>\omega_{c}$, the branch points can be integrated over without leaving
the real axis, due to Wood anomalies being known to be squareroot
singularities [29, 30, 63, 64]. They may then be dealt with using special
Gaussian quadrature, e.g. [65]. With the emergence of the trapped mode poles
at $\omega\leq\omega_{c}$, options include singularity extraction [19], and
contour deformation [7, 66]. We opt for the latter, deforming the integral
path away from the real axis and the branch cuts as shown by the solid black
line in fig. 8. Due to no poles being crossed, it follows from Cauchy’s
theorem that the value of the integral stays constant during the deformation.

Note that in order for the solution to correspond to outgoing waves, the
integration contour has to be deformed into the upper half-plane for negative
$\kappa$ and the lower half-plane for positive $\kappa$. This idea was
introduced in [67] and termed _limiting absorption principle_ by [68]. The
direction of the branch cuts was chosen with this in mind.

There exist many possible integration contours, and while there is merit in
investigating which one is most numerically efficient, we defer this to future
work ([69] discusses a systematic way of comparing contours, albeit in a
different context). Our contour of choice in fig. 8 is a sinusoid,
parameterized only by its amplitude $A$. To carry out the integration, we use
the periodic trapezoidal rule with $P_{\mathrm{asm}}$ nodes. To get an idea of
what accuracy a given $A$ and $P_{\mathrm{asm}}$ yields, we use 50 as a test
problem and investigate how accurately $\Phi$ can be reconstructed from
${\Phi_{\mathrm{p}}}_{,\kappa}$ in fig. 9. The top part of the figure then
shows the real part of total solution computed to $7$ digits of accuracy at
$\omega=2.4$, with a point source at $\mathbf{x}=(-0.2,0.1)$. The bottom part
of fig. 9 shows that the integral is sensitive to the exact value of the
amplitude, but a value between $10^{-1}$ and $1$ yields quick convergence.
With the amplitude appropriately chosen, $P_{\mathrm{asm}}\approx 100$
quadrature nodes are enough to get an answer accurate to machine precision.
This test problem, however, lacks poles in the integrand, which may put
different limits on $A$.

(a)

(b)

Figure 9: Top: Total acoustic pressure field from a point source at
$\mathbf{x}_{0}=(-0.2,0.1)$, radiating with frequency $\omega=2.4$, accurate
to $7$ digits. Bottom: convergence test of the reconstruction of
$\Phi(\mathbf{x}_{\mathrm{target}})$, at
$\mathbf{x}_{\mathrm{target}}=(0.13,0.03)$, with the array scanning method.
The sinusoidal contour shown in fig. 8 was taken, with various amplitudes $A$,
and discretized with $P_{\mathrm{asm}}$ trapezoidal quadrature nodes.

So far we only considered computation of the solution inside the central unit
cell (and above it), $\tfrac{d}{2}\leq x_{1}\leq\tfrac{d}{2}$. Using the
quasiperiodicity property, the array scanning integral 51 for a target
position in the $n$th unit cell becomes

$u^{n}(\mathbf{x})\mathrel{\mathop{\mathchar
58\relax}}=u(x_{1}+nd,x_{2})=\frac{1}{2\pi}\int_{-\pi}^{\pi}u_{\kappa}(\mathbf{x})e^{ind\kappa}\mathrm{d}\kappa.$
(52)

Figure 10 shows the array scanning integrand at a target $n=5$ unit cells away
from the center, to be compared to fig. 8. The exponential factor causes
oscillations in the real and exponential growth in the imaginary direction in
the lower (upper) half-plane for positive (negative) $n$. The exponential
growth limits the amount of contour deformation away from the real axis, and
both the growth and oscillations demand more quadrature points to maintain
constant accuracy as $n$ increases. In the limit of $n\to\infty$, however, the
solution may be computed easily as follows. Overdrawn in orange and blue for
positive and negative $n$ are the contours with which we use to make a closed
loop with the sinusoid. Taking the limit $n\to\infty$, let the orange contour
extend to infinity in the upper half-plane, and let it approach the branch cut
such that the radius of the “keyhole” around the branch point approaches zero.
In this limit, all contributions of the contour in the upper half-plane
vanish, and by Cauchy’s theorem, the value of the array scanning integral is

$\displaystyle u^{\infty}(\mathbf{x})$
$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\lim_{n\to\infty}u(x_{1}+nd,x_{2})=\lim_{n\to\infty}\frac{1}{2\pi}\int_{-\pi}^{\pi}u_{\kappa}(x_{1},x_{2})e^{in\kappa
d}\mathrm{d}\kappa$ $\displaystyle=\pm
i\Res_{\kappa=\kappa_{\mathrm{tr}}}u_{\kappa,s}(x_{1},x_{2})\lim_{n\to\pm\infty}e^{ind\kappa_{\mathrm{tr}}},$
(53)

where $\kappa_{\mathrm{tr}}$ is the position of the positive-$\kappa$ pole,
corresponding to a trapped mode.

The main claim of the previous section is a direct consequence of the above
result: at $\omega<\omega_{c}$ the only contributor to the field infinitely
far away from the source is the trapped mode at the given $\omega$. At
$\omega>\omega_{c}$, where no trapped modes exist, the field in this limit is
zero by Cauchy’s theorem, since the contour in the complex $\kappa$ plane
contains no poles. As we have seen in section 3, trapped modes at different
frequencies different vertical decay lengths and propagate at different
speeds. It is therefore of interest to compute how much of the power injected
into the system is transported away in trapped modes, as opposed to radiated
away, as a function of frequency – this is the subject of the next section.

Figure 10: Real part of the array scanning integrand 52 with the target point
$n=5$ unit cells away from the source, to be compared with fig. 8. Identical
branch cuts, poles, and contour deformation are shown in black. In orange and
blue we plot the contours we use to derive the field in the limit of
$n\to\pm\infty$, respectively.

##  5 Asymptotics of power balance

In this section we find the total power injected into the system by a single
point source, and compute the fractional power carried in trapped modes as a
function of frequency.

The total power radiated by a single point source can be derived by taking the
total flux 46 exiting a circle of radius $r$ (denoted $\Gamma_{r}$) centered
on the source, and then taking the limit $r\to 0$:

$P_{\mathrm{tot}}=\lim_{R\to
0}\imaginary\left(\int_{\Gamma_{r}}(\bar{u}_{s}+\bar{u}_{i})\partial_{n}(u_{s}+u_{i})\mathrm{d}s_{\Gamma_{r}}\right).$
(54)

Out of the four terms, only two contribute: since $u_{s}$ and its derivative
is finite everywhere, the limit of $u_{s}\partial_{n}u_{s}$ vanishes; and
$u_{i}$ only grows like $\ln r$ as $r\to 0$, therefore its contribution is
$\lim_{r\to 0}r\ln r=0$. Using the asymptotic form of $H_{0}^{(1)}(\omega r)$
as $r\to 0$ [54, Ch. 10.7], it can be shown that the remaining terms add up to

$P_{\mathrm{tot}}=\frac{1}{4}+\imaginary(u_{s}(x_{0})).$ (55)

As shown in section 4, the only contributors to the field infinitely far away
from the source along the surface are trapped modes. Therefore the power
carried in trapped modes will be calculated by considering the $n\to\pm\infty$
limit. Let us use 53 to reconstruct $u^{\infty}$ along a vertical line segment
$\Gamma_{\infty}$ at $x_{1}=\pm\tfrac{d}{2}$ (the trough of the stair),
extending from $x_{2}^{0}\mathrel{\mathop{\mathchar 58\relax}}=\tfrac{d}{2}$
to a value $x_{2}^{1}$, by which point the amplitude of the given trapped mode
has sufficiently decayed. Based on the radiation condition, we set $x_{2}^{1}$
to

$x_{2}^{1}(\kappa)\mathrel{\mathop{\mathchar
58\relax}}=-\frac{\log(10^{-4})}{\sqrt{\kappa^{2}-\omega^{2}}}.$ (56)

Then the power in trapped modes is

$P_{\mathrm{tr}}=\imaginary\left(\int_{x_{2}^{0}}^{x_{2}^{1}}\bar{u}^{\infty}\partial_{x_{1}}u^{\infty}\mathrm{d}x_{2}\right),$
(57)

which can be computed straightforwardly (note that the arbitrary phase in 53
is canceled in the above expression). Although the $x_{1}$-position of the
line segment across which the flux is calculated does not matter, numerically
the choice of the corner at $x_{1}=\pm\tfrac{d}{2}$ is ideal. This can be seen
by expanding the field in terms of Bessel functions around this point and
imposing Neumann boundary conditions; this shows that the expansion is
regular, whereas it would be singular around $x_{1}=0$. The expansion being
regular signals that the integral 57 may be carried out using Gauss–Legendre
quadrature, and no special rules are required as long as the nodes are not
within a distance comparable to the width of the smallest quadrature panel
from the boundary.

Figure 11 shows the power balance as a function of $\kappa$ on the left and
$\omega$ on the right. Since the fluxes calculated here are unitless, we show
the flux in trapped modes, the total flux, and the fraction of the flux
carried in trapped modes all in one plot. As expected from their large
vertical decay length, the least trapped modes at $\omega$, $\kappa\to 0$
carry the least flux. The highest-frequency trapped modes carry a fraction of
the total power that approaches $1$, and the increase of this fraction with
frequency or Bloch wavenumber is close to linear until an abrupt drop to zero
at $\omega>\omega_{c}$. The “critically trapped” acoustic modes are the most
efficient at capturing and transporting the input power.

(a)  (b)

Figure 11: Left: Total flux (unitless) in the system, flux in trapped modes,
and fraction of flux in trapped modes, as a function of Bloch wavenumber
$\kappa$. Right: The same quantities plotted against the frequency $\omega$ of
the Helmholtz equation.

##  6 Conclusions and future work

We described in detail a boundary integral equation method for solving the
$2$-dimensional, constant-coefficient Helmholtz equation on an exterior domain
outside of an infinite, periodic boundary with corners. We did this in the
context of acoustic scattering, with a particular emphasis on trapped modes,
which are eigenmodes of the Fredholm determinant corresponding to the
eigenvalue $0$. Using corner-refined panel quadrature on the boundary, we
built a numerical solver for an quasiperiodic array of point sources. We used
root finding to find the zeros of the Fredholm determinant and found the
dispersion relation and group velocities of all possible trapped modes.
Coupled to a simple ray model, the properties of the trapped modes were
sufficient to reconstruct the “chirp” phenomenon observed at El Castillo at
Chichen Itza, and other similar acoustics recorded at step-temples. By
integrating over the quasiperiodicity parameter, a method known as array
scanning, we computed the scattering solution from a single point source, and
extracted the limit of the acoustic pressure field infinitely far away from
the source. We showed that it is only the trapped acoustic modes that
contribute to the field at this distance. Finally, we exploited this fact to
study the fraction of injected acoustic power that ends up in trapped modes,
as a function of frequency.

With the presented numerical method in hand, several questions remain for
staircase acoustics: are there source locations close to the surface which
excite trapped modes less, or excite the left- and right-going modes
asymmetrically? Similarly, it would be of interest to extend the analysis to a
wider range of periodic boundaries that possess asymmetry.

Although dense direct linear algebra as presented here was adequate for the
studied geometry and accuracies, more complex geometries requiring more than
$10^{4}$ discretization nodes would benefit from the use of iterative solution
with FMM acceleration. Similarly, if accuracy close to machine precision was
required, more care needs to be taken when $\kappa\to\pm\pi$ and $\kappa\to 0$
due to the coalescence of poles. Here the shape of the most efficient contour
in the complex $\kappa$ plane, and its nonuniform discretization, needs
further study.

The extension to 3D doubly-periodic structures to high-order accuracy poses an
interesting challenge, because there the set of Wood-anomaly on-surface
wavevectors form curves.

## Acknowledgments

We are indebted to Eric J. Heller for suggesting that the raindrop effect is
due to the dispersion relation of trapped evanescent waves. We also thank
Charlie Epstein, Manas Rachh, Leslie Greengard, and Simon Chandler-Wilde for
their helpful comments and input. The Flatiron Institute is a division of the
Simons Foundation.

## References

  * [1] K. Busch, G. Von Freymann, S. Linden, S. Mingaleev, L. Tkeshelashvili, M. Wegener, Periodic nanostructures for photonics, Physics Reports 444 (3-6) (2007) 101–202.
  * [2] J. D. Joannopoulos, S. G. Johnson, R. D. Meade, J. N. Winn, Photonic Crystals: Molding the Flow of Light, 2nd Edition, Princeton Univ. Press, Princeton, NJ, 2008\.
  * [3] Y. Chen, Y. Zhang, A. F. Koenderink, General point dipole theory for periodic metasurfaces: magnetoelectric scattering lattices coupled to planar photonic structures, Opt. Express 25 (18) (2017) 21358–21378.
  * [4] G. Ji, J. Huber, Recent progress in acoustic metamaterials and active piezoelectric acoustic metamaterials-a review, Applied Materials Today 26 (2022) 101260.
  * [5] J. Lin, H. Zhang, An integral equation method for numerical computation of scattering resonances in a narrow metallic slit, J. Comput. Phys. 385 (2019) 75–105.
  * [6] B. Munk, G. Burrell, 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 transactions on antennas and propagation 27 (3) (1979) 331–343.
  * [7] Q.-Q. He, B.-Z. Wang, Radiation patterns synthesis for a conformal dipole antenna array, Progress In Electromagnetics Research 76 (2007) 327–340.
  * [8] I. Herrero-Durá, A. Cebrecos, R. Picó, V. Romero-García, L. M. García-Raffi, V. J. Sánchez-Morcillo, Sound absorption and diffusion by 2d arrays of helmholtz resonators, Applied Sciences 10 (5) (2020) 1690.
  * [9] N. F. Declercq, C. S. Dekeyser, Acoustic diffraction effects at the Hellenistic amphitheater of Epidaurus: Seat rows responsible for the marvelous acoustics, The Journal of the Acoustical Society of America 121 (4) (2007) 2011–2022.
  * [10] S. Shipman, Resonant scattering by open periodic waveguides, Vol. 1 of Progress in Computational Physics (PiCP), Bentham Science Publishers, 2010, Ch. 2, pp. 7–50.  doi:10.2174/978160805150211001010007.
  * [11] S. N. Chandler-Wilde, I. G. Graham, S. Langdon, E. A. Spence, Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering, Acta numerica 21 (2012) 89–305.
  * [12] E. L. Richards, H. Song, W. Hodgkiss, Acoustic scattering comparison of kirchhoff approximation to rayleigh-fourier method for sinusoidal surface waves at low grazing angles, The Journal of the Acoustical Society of America 144 (3) (2018) 1269–1278.
  * [13] N. Tsingos, S. Lefebvre, C. Dachsbacher, M. Dellepiane, Extending geometrical acoustics to highly detailed architectural environments, in: 19th International Congress on Acoustics (ICA 2007), 2007, p. 6.
  * [14] J. A. C. Calleja, N. F. Declercq, The acoustic raindrop effect at Mexican pyramids: The architects’ homage to the rain god Chac?, Acustica 95 (2009) 849–856.
  * [15] D. Colton, R. Kress, Inverse acoustic and electromagnetic scattering theory, 2nd Edition, Vol. 93 of Applied Mathematical Sciences, Springer-Verlag, Berlin, 1998.
  * [16] M. S. Howe, Acoustics of fluid-structure interactions, Cambridge University Press, 1998.
  * [17] D. Colton, R. Kress, Integral Equation Methods in Scattering Theory, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1983.  doi:10.1137/1.9781611973167.
  * [18] M. Kaltenbacher, Computational acoustics, Springer, 2018.
  * [19] I. Rana, N. Alexopoulos, Current distribution and input impedance of printed dipoles, IEEE Transactions on Antennas and Propagation 29 (1) (1981) 99–105.
  * [20] F. Capolino, D. R. Jackson, D. R. Wilton, Mode excitation from sources in two-dimensional EBG waveguides using the array scanning method, IEEE Microwave and Wireless Components Letters 15 (2) (2005) 49–51.
  * [21] A. Lechleiter, R. Zhang, A convergent numerical scheme for scattering of aperiodic waves from periodic surfaces based on the Floquet–Bloch transform, SIAM J. Numer. Anal. 55 (2) (2017) 713–736.
  * [22] R. Zhang, Numerical methods for scattering problems in periodic waveguides, Numer. Math. 148 (2021) 959–996.  doi:10.1007/s00211-021-01229-0.
  * [23] A.-S. Bonnet-BenDhia, F. Starling, Guided waves by electromagnetic gratings and non-uniqueness examples for the diffraction problem, Math. Meth. Appl. Sci. 17 (1994) 305–338.
  * [24] C. H. Wilcox, Scattering theory for diffraction gratings, Applied Mathematical Sciences, volume 46, Springer-Verlag, 1984.
  * [25] D. V. Evans, C. M. Linton, Edge waves along periodic coastlines, Q. J. Mech. Appl. Math. 46 (1993) 643–656.
  * [26] V. Y. Gotlib, Solutions of the Helmholtz equation, concentrated near a plane periodic boundary, J. Math. Sci. 102 (2000) 4188–4194.  doi:10.1007/BF02673850.
  * [27] C. L. Epstein, Solving the transmission problem for open wave-guides. II outgoing estimates (2023).  arXiv:2310,05816.
  * [28] R. W. Wood, On a remarkable case of uneven distribution of light in a diffraction grating spectrum, Philos. Mag. 4 (1902) 396–408.
  * [29] U. Fano, The theory of anomalous diffraction gratings and of quasi-stationary waves on metallic surfaces (Sommerfeld’s waves), J. Opt. Soc. Am. 31 (3) (1941) 213–222.
  * [30] A. Hessel, A. Oliner, A new theory of Wood’s anomalies on optical gratings, Appl. Opt. 4 (10) (1965) 1275–1297.
  * [31] T. Arens, S. N. Chandler-Wilde, J. A. DeSanto, On integral equation and least squares methods for scattering by diffraction gratings, Commun. Comput. Phys. 1 (2006) 1010–42.
  * [32] O. P. Bruno, M. C. Haslam, Efficient high-order evaluation of scattering by periodic surfaces: deep gratings, high frequencies, and glancing incidences, J. Opt. Soc. Am. A 26 (3) (2009) 658–668.
  * [33] A. Barnett, L. Greengard, A new integral representation for quasi-periodic scattering problems in two dimensions, BIT Numerical mathematics 51 (2011) 67–90.
  * [34] O. P. Bruno, B. Delourme, Rapidly convergent two-dimensional quasi-periodic Green function throughout the spectrum—including Wood anomalies, J. Comput. Phys. 262 (2014) 262–290.
  * [35] M. H. Cho, A. H. Barnett, Robust fast direct integral equation solver for quasi-periodic scattering problems with a large number of layers, Optics express 23 (2) (2015) 1775–1799.
  * [36] R. Kress, V. Maz’ya, V. Kozlov, Linear integral equations, Vol. 82, Springer, 1989\.  doi:10.1007/978-1-4612-0559-3.
  * [37] R. Zhang, A high order numerical method for scattering from locally perturbed periodic surfaces, SIAM Journal on Scientific Computing 40 (4) (2018) A2286–A2314.
  * [38] G. Fairweather, A. Karageorghis, The method of fundamental solutions for elliptic boundary value problems, Advances in Computational Mathematics 9 (1998) 69–95.
  * [39] A. H. Cheng, Y. Hong, An overview of the method of fundamental solutions—solvability, uniqueness, convergence, and stability, Engineering Analysis with Boundary Elements 120 (2020) 118–152.
  * [40] A. H. Barnett, T. Betcke, Stability and convergence of the method of fundamental solutions for helmholtz problems on analytic domains, Journal of Computational Physics 227 (14) (2008) 7003–7026.
  * [41] C. J. Alves, S. S. Valtchev, Numerical comparison of two meshfree methods for acoustic wave scattering, Engineering Analysis with Boundary Elements 29 (4) (2005) 371–382.
  * [42] R. Petit (Ed.), Electromagnetic Theory of Gratings, Vol. 22 of Topics in Current Physics, Springer-Verlag,Heidelberg, 1980.
  * [43] R. Millar, The rayleigh hypothesis and a related least-squares solution to scattering problems for periodic surfaces and other scatterers, Radio Science 8 (8-9) (1973) 785–796.
  * [44] W. C. Meecham, On the use of the kirchhoff approximation for the solution of reflection problems, Journal of Rational Mechanics and Analysis 5 (2) (1956) 323–334.
  * [45] C. M. Linton, The Green’s function for the two-dimensional Helmholtz equation in periodic domains, J. Eng. Math. 33 (1998) 377–402.
  * [46] C. M. Linton, Lattice sums for the Helmholtz equation, SIAM review 52 (4) (2010) 630–674.
  * [47] K. Yasumoto, K. Yoshitomi, Efficient calculation of lattice sums for free-space periodic Green’s function, IEEE Transactions on Antennas and Propagation 47 (6) (1999) 1050–1055.
  * [48] J. L. Tsalamengas, Gauss–jacobi quadratures for weakly, strongly, hyper-and nearly-singular integrals in boundary integral equation methods for domains with sharp edges and corners, Journal of Computational Physics 325 (2016) 338–357.
  * [49] T. A. Driscoll, L. N. Trefethen, Schwarz-Christoffel mapping, Vol. 8, Cambridge University Press, 2002.
  * [50] J. Bremer, V. Rokhlin, I. Sammis, Universal quadratures for boundary integral equations on two-dimensional domains with corners, Journal of Computational Physics 229 (22) (2010) 8259–8280.
  * [51] J. Helsing, R. Ojala, Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning, Journal of Computational Physics 227 (20) (2008) 8820–8840.
  * [52] A. Gopal, L. N. Trefethen, New Laplace and Helmholtz solvers, Proceedings of the National Academy of Sciences 116 (21) (2019) 10223–10225.
  * [53] L. Zhao, A. Barnett, Robust and efficient solution of the drum problem via nystrom approximation of the fredholm determinant, SIAM Journal on Numerical Analysis 53 (4) (2015) 1984–2007.
  * [54] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark (Eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, 2010.
URL http://dlmf.nist.gov

  * [55] L. N. Trefethen, Approximation Theory and Approximation Practice, Extended Edition, SIAM, 2019.
  * [56] L. Greengard, K. L. Ho, J.-Y. Lee, A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions, Journal of Computational Physics 258 (2014) 738–751.
  * [57] K. E. Atkinson, The numerical solution of integral equations of the second kind, Cambridge University Press, 1997.
  * [58] W. H. Press, Numerical recipes 3rd edition: The art of scientific computing, Cambridge university press, 2007.
  * [59] A. Kirsch, Uniqueness theorems in inverse scattering theory for periodic structures, Inv. Probs. 10 (1) (1994) 145–152.  doi:10.1088/0266-5611/10/1/011.
  * [60] L. Zhao, A. H. Barnett, Robust and efficient solution of the drum problem via Nyström approximation of the Fredholm determinant, SIAM J. Numer. Anal. 53 (4) (2015) 1984–2007.
  * [61] O. Bühler, A Brief Introduction to Classical, Statistical, and Quantum Mechanics, Courant Lecture Notes, AMS, 2006.
  * [62] N. F. Declercq, J. Degrieck, R. Briers, O. Leroy, A theoretical study of special acoustic effects caused by the staircase of the el castillo pyramid at the maya ruins of chichen-itza in mexico, The Journal of the Acoustical Society of America 116 (6) (2004) 3328–3335.
  * [63] B. Bolotovskii, A. Lebedev, On threshold phenomena in classical electrodynamics, Sov. Phys. JETP 26 (4) (1968) 784–786.
  * [64] C. C. Wojcik, H. Wang, M. Orenstein, S. Fan, Universal behavior of the scattering matrix near thresholds in photonics, Physical Review Letters 127 (27) (2021) 277401.
  * [65] D. Huybrechs, R. Cools, On generalized Gaussian quadrature rules for singular and nearly singular integrals, SIAM journal on numerical analysis 47 (1) (2009) 719–739.
  * [66] G. Lovat, R. Araneo, S. Celozzi, Dipole excitation of periodic metallic structures, IEEE transactions on antennas and propagation 59 (6) (2011) 2178–2187.
  * [67] W. v. Ignatowsky, Reflexion elektromagnetisches Wellen an einem Draht, Annalen der Physik 323 (13) (1905) 495–522.
  * [68] A. Sveshnikov, The radiation principle, in: Dokl. Akad. Nauk SSSR, Vol. 73, 1950, pp. 917–920.
  * [69] A. Hochman, Y. Leviatan, A numerical methodology for efficient evaluation of 2D Sommerfeld integrals in the dielectric half-space problem, IEEE transactions on antennas and propagation 58 (2) (2009) 413–431.