We investigate eigenvalue problems for the planar Helmholtz equation in open systems with a high order of rotational symmetry. The resulting solutions have similarities with the whispering gallery modes exploited in photonic micro-resonators and elsewhere, but unlike these do not necessarily require a surrounding material boundary, with confinement instead resulting from the geometry of a series of inclusions arranged in a ring. The corresponding fields exhibit angular quasi-periodicity reminiscent of Bloch waves, and hence we refer to them as whispering Bloch modes (WBMs). We show that if the geometry of the system is slightly perturbed such that the rotational symmetry is broken, modes with asymmetric field patterns can be observed, resulting in field enhancement and other potentially desirable effects. We investigate the WBMs of two specific geometries first using expansion methods and then by applying a two-scale asymptotic scheme.

1. Introduction

The whispering gallery phenomenon can famously be observed in the auditorium of St. Paul’s Cathedral, where a message whispered next to the surrounding wall of the circular viewing gallery can be heard clearly at any other point on the circumference, while remaining unintelligible to an observer sitting away from the wall. A mathematical explanation was provided by Lord Rayleigh in 1910 [1] and since then analogous effects, ubiquitous in wave systems with circular or spherical cavities, have been well studied [2]. In a modern setting, whispering gallery modes (WGMs) in electromagnetic wave systems have found application in spectroscopy, microdisk lasers [3], biosensors [4] and experimental testing of nonlinear optics [5]. In most instances, the excitation is partially confined due to a material boundary, and as a result the solutions are in fact quasi-modes; assuming harmonic time dependence $\exp (- i ω t)$ is understood, they are characterized by complex frequencies with Im(ω)<0, and hence decay in time as energy is radiated to infinity. A figure of merit is the so-called radiative quality factor Q=Re(ω)/|2Im(ω)|.

In this paper, we consider open systems with high orders of rotational symmetry. It is straightforward to prove that the quasi-modes of such systems exhibit discrete Bloch-like angular periodicity and hence we refer to them as WBMs. Unlike WGMs, these quasi-modes do not require a surrounding material boundary to achieve confinement and this can instead be induced by a geometric structure arranged periodically around a circle. To clearly distinguish WBMs from WGMs, in figure 1, we show an illustrative example displaying four apparently very similar modes. The figure depicts TE-polarized eigenstates of a high-index cylinder in vacuum, in which we place an array of spokes with Neumann boundary conditions imposed upon them (representing perfect conductors in this polarization). Irrespective of the presence of the spokes, a WGM is found due to the material mismatch between the cylinder and background, as shown in figure 1a,b. In the case with spokes, however, there is an additional mode, shown in figure 1c, that is created by the periodicity of the geometry; this second mode persists even if the materials are identical, as shown in figure 1d. We refer to the latter as a WBM. Even without the confining effect of the material boundary, in this case, the Q-factor is several orders of magnitude greater for the WBM than for the WGM.

Figure 1. Whispering gallery modes versus whispering Bloch modes. Frame (a) shows a conventional WGM for a high refractive index (n_r=2) cylinder in a low refractive index (n_r=1) background. In (b), we demonstrate that the same solution exists if the cylinder is embedded with N=30 radial Neumann spokes. Frame (c) shows a WBM that exists in the latter system and is a result of the spokes themselves; a similar mode exists even if the material boundary is removed, as shown in frame (d). In each case, the real part of the field is plotted, normalized to the colour bar shown. Note that the theoretical Q-factors are far higher for the lower panels, suggesting a greater degree of confinement.

In the remainder of this article, we consider cyclic arrays of inclusions embedded within a homogeneous background and therefore WGMs cannot exist. We study systems governed by the planar Helmholtz equation with Neumann boundary conditions and use two geometries to illustrate WBMs, each of which can be treated analytically using an appropriate expansion method. The first example, an arrangement of equally spaced radial spokes, is instructive as the solutions bear close resemblance to conventional WGMs as shown in figure 1. This eigenvalue problem is solved using matched eigenfunction expansions. The second example is an array of isolated circular inclusions, which is closely related to a linear array of cylinders that is well known to support Rayleigh–Bloch waves [6], but now bent into a closed curve. In this case, the methodology of multipoles is used instead. A detailed study of resonances in rings of circular inclusions was undertaken in [7], and our work is complementary to this, extending the regime of interest to larger numbers of inclusions, for which new and interesting effects are observed. More recently, similar systems have been studied as models for acoustic and electric Faraday cages; the article by Chapman et al. [8] focuses on shielding effects in the static and quasi-static regime, while Martin [9] develops a wire model using Foldy’s method.

Much of the motivation for this paper comes from the extensive literature on Rayleigh–Bloch surface waves on linear or planar arrays; solutions for the Helmholtz equation have been studied for comb-like structures [10,11], and also for arrays of rectangular inclusions [12], cylinders [6] and more irregular shapes [13]. General existence results have been derived for surfaces with Neumann boundary conditions [14], and for this reason, we focus our attention on these rather than the alternative Dirichlet conditions for which equivalent modes do not exist [15]. In the context of electromagnetic waves, guided modes are observed in linear or planar arrays of dielectric spheres [16] and also on conducting surfaces with periodic corrugations that support spoof surface plasmon polaritons [17,18]. It was recently shown that the quasi-modes of a grooved conducting disc give rise to magnetic localized plasmons [19]. Furthermore, there is an older literature on electromagnetic gratings showing surface waves guided by corrugations [20,21].

At a discrete set of frequencies, the angular periodicity of a system’s WBMs will match that of the rotationally symmetric structure or will be a half-integer multiple thereof. In these cases, the solutions satisfy periodic or antiperiodic boundary conditions across a single wedge-shaped repeating cell, and these are the standing wave eigenmodes of the system. Analogous standing wave solutions are found in structures with Bravais lattice periodicity [22], and this suggests that an asymptotic method based on a cyclic analogue of high-frequency homogenization [23] would be well suited to studying the systems in question. We investigate this and find, at complex frequencies close to those of the standing wave solutions, quasi-modes exist with wide-angle modulation, and we construct asymptotic expressions for the field patterns and associated complex frequencies; these modulation effects are not observed in WGM systems. We then go on to investigate the effect of slightly perturbing the geometry in a WBM system, and found that the asymptotic method can be adapted to account for this, predicting the perturbed field patterns and frequencies.

Further to the methods mentioned above, we employ finite-element method (FEM) modelling, first as a numerical check, and later in conjunction with the asymptotic theory we develop. With this approach we necessarily truncate the domain to a finite size, and particularly in the case of leakier modes with low Q-factors this can introduce a significant error in the calculated eigenvalues. In terms of solving the eigenvalue problem, the analytical methods are preferable as the solutions are written naturally in bases that satisfy the outgoing wave condition, and hence no such truncation is required.

We begin in §2 with the general formulation of the problem, along with a proof of the angular quasi-periodicity condition satisfied by the solutions. We then turn to the two examples that we primarily consider: the array of radial spokes and the array of isolated holes. In both cases, we use analytical methods to generate systems of equations that we solve numerically to generate discrete dispersion plots. An asymptotic method is advanced in §3 that is insightful for analysing WBMs with wide-angle modulation that occur at frequencies close to the standing wave eigenfrequencies of the system. In §4, we adapt this asymptotic method to analyse a case where the geometry of the system is slightly perturbed, leading to asymmetric field effects. Comparisons with FEM simulation are made to demonstrate the accuracy and simplicity of the approach. Finally, we draw together some concluding remarks in §5.

2. Prototype system

We pose an eigenvalue problem for the planar Helmholtz equation with Neumann inclusions, seeking $u : R^{2} ∖ \cup C_{j} \to C$ , where $u \in C^{2} (R^{2} ∖ \cup {\bar{C}}_{j})$ and $Ω \in C$ satisfy

\begin{aligned} (Δ + Ω^{2}) u (x) = 0, x \in R^{2} ∖ \cup {\bar{C}}_{j}, \\ n \cdot \nabla u (x) = 0, x \in \cup \partial C_{j} \\ and & \frac{\partial u}{\partial r} \to i Ω u, r \to \infty . \end{aligned}}

2.1

Here

C_{j}

for

j = 1, \dots, \tilde{N}

refers to a set of bounded inclusions arranged to have rotational symmetry of order N. Note that

\tilde{N}

may be greater or less than N, and a simple case to consider would be that of

\tilde{N} = N

identical inclusions arranged periodically around a ring. The unit normal vector at an inclusion boundary is denoted by n, and r=|x|. In particular, we consider the two geometries shown in figure 2 so that concrete solutions can be found.

Figure 2. Geometries under consideration; (a) a system of N equally spaced radial spokes of unit length and infinitesimal width and (b) a ring of N circular holes of radius a, whose centres are equally spaced around a ring of unit radius; we take $a < \sin (π / N)$ so that the holes do not overlap.

Before solving the eigenvalue problem for the two geometries under consideration, we derive a quasi-periodicity condition analogous to Bloch’s theorem for Bravais lattice structures. Consider a simple eigenvalue Ω of system (2.1), along the associated field u(r,ϕ). We introduce a discrete rotation operator $\hat{R}$ such that for any suitably defined function, $\hat{R} f (r, ϕ) = f (r, ϕ + 2 π / N)$ . The rotational symmetry of (2.1) means that the field $\hat{R} u (r, ϕ)$ solves an identical eigenvalue problem to u(r,ϕ), and since the eigenvalue is simple, the two must be identical up to multiplication by a complex number, i.e. $\hat{R} u (r, ϕ) = α u (r, ϕ)$ for $α \in C$ . If we perform the same transformation N times, the cyclic continuity of the problem requires that $α^{N} u (r, ϕ) = {\hat{R}}^{N} u (r, ϕ) = u (r, ϕ)$ and hence α=e^2mπi/N is one of the Nth roots of unity. The rotational analogue of Bloch’s theorem is thus

u (r, ϕ + \frac{2 j π}{N}) = e^{2 j m π i / N} u (r, ϕ), \forall j, m \in Z .

2.2

We do not consider the case of degenerate eigenvalues in this paper, and this possibility is left as an open question. The above result leads to a discrete Bloch-like spectrum in complex frequency space for m=0,1,…,N−1. In the examples chosen, symmetry under the transformation ϕ→−ϕ allows us to consider a reduced domain m=0,1,…,⌊N/2⌋, analogous to the irreducible Brillouin zone in conventional Bloch theory [24]. The spectra for both systems are calculated using expansion methods. We begin with the geometry of figure 2a.

(a) Neumann spokes: matched eigenfunction expansion method

We consider the infinite wedge-shaped repeating cell as in figure 2a. Supposing that one of the spokes lies on the half-line ϕ=0, we expand the field in the two regions shown in figure 3 using appropriate Fourier–Bessel bases. In particular, the Neumann condition must be satisfied by each term in $D_{1}$ , and the quasi-periodicity condition must be satisfied by each term in $D_{2}$ . The expansions are thus given by

u (r, ϕ) = {\begin{cases} \sum_{n = 0}^{\infty} {\hat{A}}_{n} J_{n N / 2} (Ω r) \cos (\frac{n N ϕ}{2}), & (r, ϕ) \in D_{1}, \\ \sum_{n = - \infty}^{\infty} {\hat{B}}_{n} H_{m + n N} (Ω r) e^{i (m + n N) ϕ}, & (r, ϕ) \in D_{2}, \end{cases}

2.3

where J and H denote the Bessel and the Hankel functions of the first kind, respectively. To ensure that the field has the required smoothness, the expressions for u and its derivative ∂u/∂r from these two expansions are equated at r=1. We multiply the resulting equations by

\cos (p N ϕ / 2)

for

p \in Z

and integrate over the azimuthal variable ϕ. Using the orthogonality of the cosine functions, along with the relation

\begin{aligned} β_{p n} & \equiv \int_{0}^{2 π / N} e^{i (m + n N) ϕ} \cos (\frac{p N ϕ}{2}) d ϕ \\ = {\begin{cases} \frac{π}{N}, & m + n N = \pm \frac{p N}{2}, \\ \frac{i (m + n N) (1 - e^{i m 2 π / N} {(- 1)}^{p})}{(m + n N + p N / 2) (m + n N - p N / 2)}, & otherwise, \end{cases} \end{aligned}

2.4

we are able to eliminate the coefficients

{\hat{A}}_{n}

, leaving an infinite system of equations:

0 = \sum_{n = - \infty}^{\infty} {\hat{P}}_{p n} {\hat{B}}_{n} \forall p \in Z,

2.5

where

{\hat{P}}_{p n} = μ_{p n}^{(1)} - μ_{p n}^{(2)}

and

μ_{p n}^{(1)} = [J_{p N / 2 - 1} (Ω) - J_{p N / 2 + 1} (Ω)] H_{m + n N} (Ω) β_{p n}

2.6

and

μ_{p n}^{(2)} = [H_{m + n N - 1} (Ω) - H_{m + n N + 1} (Ω)] J_{p N / 2} (Ω) β_{p n} .

2.7

Figure 3. Geometry of the cell $S_{0}$ . Separate expansions are used in the regions $D_{1} = {(r, ϕ) : r \in (0, 1], ϕ \in [- π / N, π / N]}$ and $D_{2} = {(r, ϕ) : r \in [1, \infty), ϕ \in [- π / N, π / N]}$ .

We note that in the study of Rayleigh–Bloch surface waves on a comb-like grating, a particularly elegant method involving residue calculus can be employed to solve the system of equations analogous to (2.5) [11,20], but unfortunately the properties of the Bessel functions in the present case prevent the same approach from being used straightforwardly here. The residue calculus method makes explicit use of the expected square-root singularity of the solution at the tip of each spoke, and a result that we can borrow readily is that the same singularity implies that ${\hat{B}}_{n} \sim O (n^{- 3 / 2})$ as $n \to \infty$ . Rather than pursuing an analytical method, we truncate the infinite system (2.5) to a finite system of 2M+1 equations with −M≤n≤M, 0≤p≤2M, which we then solve numerically. A problem arises here as the higher order Bessel and Hankel functions have very small and very large magnitudes, respectively, so the resulting matrix is ill-conditioned. To remedy this, we introduce matrix elements given by $P_{p n} = {\hat{P}}_{p n} / (J_{p N / 2} (Ω) H_{m + n N} (Ω))$ , which satisfy $0 = \sum_{n} P_{p n} B_{n}$ $\forall p \in Z$ for the same values of $Ω \in C$ as before, but do not blow up as |n| or |p| become large. This resolves the ill-conditioning problem, but the magnitudes of the Bessel and Hankel functions themselves are still prohibitive for relatively modest values of M, resulting in floating-point underflow/overflow in some of the matrix elements, and failure of the scheme. In these cases, we use the asymptotic forms [25]:

J_{ζ} (z) \sim \frac{1}{\sqrt{2 π ζ}} {(\frac{e z}{2 ζ})}^{ζ} and H_{ζ} (z) \sim - i \sqrt{\frac{2}{π ζ}} {(\frac{e z}{2 ζ})}^{- ζ} as ζ \to + \infty along the real axis,

2.8

along with the identity J_−ζ(z)=(−1)^ζJ_ζ(z), and similarly for H_−ζ(z), in place of the corresponding functions in P_pn, and after simplification these are used in the truncated matrix equation PB=0, which has non-trivial solutions iff

det (P) = 0

. For a given value of m, there will be a spectrum

{Ω} \subset C

corresponding to the quasi-modes of the open systems (figure 4), and the expansion coefficients are recovered from the null space vector B.

(b) Neumann holes: multipole method

Other geometries are not generally well suited to eigenfunction matching, but the prototypical circular inclusion is ideally suited to multipole techniques [26,27]. Such a method was previously applied in the context of a circular ring in [7]. To solve the eigenvalue problem for the geometry of figure 2b, we define local systems of polar coordinates centred at each of the j holes (r_j,ϕ_j), as shown in figure 5. In the vicinity of the 0th hole $C_{0}$ , the field is expanded in terms of regular and outgoing waves as

u (r_{0}, ϕ_{0}) = \sum_{μ = - \infty}^{\infty} [A_{μ}^{0} J_{μ} (Ω r_{0}) + E_{μ}^{0} H_{μ} (Ω r_{0})] e^{i μ ϕ_{0}},

2.9

where

A_{μ}^{0} J_{μ}^{'} (Ω a) + E_{μ}^{0} H_{μ}^{'} (Ω a) = 0

\forall μ \in Z

ensures that the Neumann condition is satisfied on the inclusion. This expansion is valid within an annulus surrounding

C_{0}

, extending as far as the nearest point on each of the neighbouring holes

C_{1}

and

C_{N - 1}

. A second way to express the field is due to Wijngaard [28], who argued that in the absence of sources, the field existing between a number of scattering bodies must be the sum of the outgoing waves from each body. In the present case, this leads to the expansion

u (r_{0}, \dots, r_{N - 1}, ϕ_{0}, \dots, ϕ_{N - 1}) = \sum_{j = 0}^{N - 1} \sum_{ν = - \infty}^{\infty} E_{ν}^{j} H_{ν} (Ω r_{j}) e^{i ν ϕ_{j}},

2.10

where the jth term in the first summation is given in terms of coordinates centred at the jth hole, and hence u is expressed as a function of 2N (dependent) variables. This expansion is valid everywhere in

R^{2} ∖ \cup C_{i}

, and automatically satisfies the outgoing radiation condition as

r \to \infty

. Quasi-periodicity implies that

E_{ν}^{j} = E_{ν}^{0} e^{2 j m π i / N}

\forall ν \in Z

. We use Graf’s addition theorem [29] to express every term in (2.10) in the coordinate system of the 0th hole:

\begin{aligned} u (r_{0}, ϕ_{0}) & = \sum_{μ = - \infty}^{\infty} E_{μ}^{0} H_{μ} (Ω r_{0}) e^{i μ ϕ_{0}} \\ - \sum_{j = 1}^{N - 1} \sum_{ν = - \infty}^{\infty} A_{ν}^{0} \frac{J_{ν}^{'} (Ω a)}{H_{ν}^{'} (Ω a)} \sum_{μ = - \infty}^{\infty} i^{ν - μ} e^{i (2 m - μ - ν) π j / N} H_{μ - ν} (Ω r_{0 j}) J_{μ} (Ω r_{0}) e^{i μ ϕ_{0}}, \end{aligned}

2.11

where r_0j is the distance between the centres of the 0th and jth holes. This expansion has the same domain of validity as (2.9), and indeed the two expressions must be equivalent. We equate the two expansions term by term, which yields,

\forall μ \in Z

\sum_{ν = - \infty}^{\infty} A_{ν}^{0} [i^{ν - μ} \frac{J_{ν}^{'} (Ω a)}{H_{ν}^{'} (Ω a)} \sum_{j = 1}^{N - 1} e^{i (2 m - μ - ν) π j / N} H_{μ - ν} (Ω r_{0 j}) + δ_{μ ν}] = 0.

2.12

For numerical calculation, the infinite system (2.12) is truncated to a finite system of 2M+1 equations with −M≤μ,ν≤M. The corresponding matrix equation PA=0 has non-trivial solutions iff

det (P) = 0

. For a given value of m, there will be a spectrum

{Ω} \subset C

corresponding to the quasi-modes of the open systems (figures 6 and 7), and the expansion coefficients are recovered from the null space vector A.

Figure 5. Geometry of the cell S_j. In order for the holes $C_{j}$ and $C_{j \pm 1}$ to be non-overlapping, we require that $a < \sin (π / N)$ .

Figure 6. (a) Contour plot of $\log | det (P) |$ in the complex frequency plane for a unit ring of N=20 Neumann holes with radius $0.8 \sin (π / 20)$ for standing wave solutions with azimuthal quantum number m=10. Note the presence of two modes very close to the real Ω axis, corresponding to confined quasi-modes. (b,c) The real part of the field u for these confined WBMs at Ω≈14.12 and Ω≈18.90, respectively, normalized to the colour bar shown.

Figure 7. (a) Dispersion of the lowest branch of WBMs for a system of N=60 holes with radius $0.8 \sin (π / 60)$ . The blue crosses (red circles) show the real (imaginary) parts of the complex frequencies calculated using the method outlined in this section, and the dashed line shows the asymptotic dispersion relation as calculated from (3.10) in §3. (b–d) The real part of the field u for modes with m=30,29,28, respectively, normalized to the colour bar shown. Note that the standing wave solution is antiperiodic (m=N/2) in this case.

3. Asymptotic analysis

The existence of standing wave WBMs, which are periodic or antiperiodic across a single wedge-shaped elementary cell, is interesting both in terms of application and analysis. In this section, we consider large values of N and apply an asymptotic perturbation scheme to these standing wave solutions that allows us to calculate the field pattern, as well as the complex frequencies, of wide-angle modulated quasi-modes of the system. Later, we shall see that the same method is particularly useful in analysing cases where the geometry of the system is slightly perturbed, potentially leading to field-enhancement and other asymmetric effects.

The asymptotic method, outlined in figure 8, is based on high-frequency homogenization [23], which was used in the context of Rayleigh–Bloch modes on linear gratings in [10,30]. The two-scale approach relies on the disparity between the small angle subtended by the wedge-shaped elementary cell, and the large angle subtended by the underlying modulation of the quasi-mode. Motivated by this, let us consider the system (2.1), but restrict our analysis to the (arbitrary) truncated cell $S_{0} = {(r, ϕ) : r \in (0, r_{+}], ϕ \in [- π / N, π / N]} ∖ \cup C_{j}$ , where r₊≫1. We impose the outgoing wave condition ∂u/∂r=(iΩ−1/2r)u at r=r₊, and for a sufficiently large value of r₊, solutions of this problem will converge to those of the infinite problem. The standing wave WBMs satisfy periodic/antiperiodic boundary conditions on opposing sides of the cell:

u |_{ϕ = π / N} = \pm u |_{ϕ = - π / N}, {\frac{\partial u}{\partial ϕ} |}_{ϕ = π / N} = \pm {\frac{\partial u}{\partial ϕ} |}_{ϕ = - π / N},

3.1

where ‘−’ is possible only if N is even. In order to exploit the scale disparity, we introduce a periodic angular variable satisfying φ∈[−π/N,π/N] in each cell, along with a slowly-varying angular variable θ=ηϕ, where η≪1. In the asymptotic limit η→0, any function of the original angular variable ϕ may be treated instead as a function of the two variables φ and θ, which are considered to be independent. The value of the parameter η is dynamic; its value is of the order of the ratio between the short-scale field variation and the long-scale Bloch modulation, but at this stage it need not be defined explicitly. The gradient and Laplacian operator are written in terms of the new variables using the chain rule as

\nabla \to \hat{r} \frac{\partial}{\partial r} + \frac{\hat{ϕ}}{r} (\frac{\partial}{\partial φ} + η \frac{\partial}{\partial θ}) \equiv \nabla_{r, φ} + \frac{\hat{ϕ}}{r} η \frac{\partial}{\partial θ}

3.2

and

Δ \to \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial}{\partial r}) + \frac{1}{r^{2}} (\frac{\partial^{2}}{\partial φ^{2}} + 2 η \frac{\partial^{2}}{\partial φ \partial θ} + η^{2} \frac{\partial^{2}}{\partial θ^{2}}) \equiv Δ_{r, φ} + \frac{1}{r^{2}} (2 η \frac{\partial^{2}}{\partial φ \partial θ} + η^{2} \frac{\partial^{2}}{\partial θ^{2}}) .

3.3

Following the approach of Craster et al. [23], we seek solutions in the form of an asymptotic ansatz:

\begin{aligned} u (r, φ, θ) & = u_{0} (r, φ, θ) + η u_{1} (r, φ, θ) + η^{2} u_{2} (r, φ, θ) + \dots \\ and Ω^{2} & = Ω_{0}^{2} + η Ω_{1}^{2} + η^{2} Ω_{2}^{2} + \dots, \end{aligned}}

3.4

noting that the field u is now considered to be a function of both angular variables that are treated as independent. The resulting hierarchy is solved with periodic/antiperiodic boundary conditions on the short scale. At leading order, we have

O (1) : {\begin{cases} (Δ_{r, φ} + Ω_{0}^{2}) u_{0} = 0, & (r, ϕ) \in S_{0}, \\ n \cdot \nabla_{r, φ} u_{0} = 0, & (r, ϕ) \in \cup \partial C_{j} \cap {\bar{S}}_{0} . \end{cases}

3.5

Figure 8. Outline of the asymptotic scheme applied to the second-lowest branch of WBMs for a ring of N=60 holes. (a) The real part of U₀ for a standing wave quasi-mode with antiperiodic boundary conditions, (b) the corresponding real part of the function U₁ and (c) the real part of the asymptotic field u₀+ηu₁ for a mode with Δm=1.

This problem is most straightforwardly solved using an FEM solver such as Comsol [31]. Assuming the eigenvalue is simple, the solution consists of a standing wave eigenfunction U₀(r,φ), multiplied by a coefficient that generally depends on the independent long-scale angular coordinate θ, so we have u₀(r,φ,θ)=f₀(θ)U₀(r,φ) for an arbitrary function f₀. At this point, we ensure that the value of r₊ is sufficiently large that the solution has converged to that of the infinite problem (2.1), which is checked by comparing values of Ω₀ with those calculated using the expansion methods in §2. Since we are interested in spatially confined quasi-modes, it is assumed that Ω₀ has only a small imaginary part. Specifically, the condition that r₊<1/|2Im (Ω₀)| should be satisfied, which ensures the hierarchy does not break down as discussed in appendix A(b).

The next order in the hierarchy poses a forced problem:

O (η) : {\begin{cases} (Δ_{r, φ} + Ω_{0}^{2}) u_{1} = - (\frac{2}{r^{2}} \frac{\partial^{2}}{\partial φ \partial θ} + Ω_{1}^{2}) u_{0}, & (r, ϕ) \in S_{0}, \\ n \cdot \nabla_{r, φ} u_{1} = - n \cdot \hat{ϕ} \frac{1}{r} \frac{\partial u_{0}}{\partial θ}, & (r, ϕ) \in \cup \partial C_{j} \cap {\bar{S}}_{0} . \end{cases}

3.6

At this point, we derive a compatibility condition that implies Ω₁=0 (appendix A(a)). With this in place, the forced system is solved, yielding a solution of the form u₁(r,φ,θ)=f₁(θ)U₀(r,φ)+f₀′(θ)U₁(r,φ), where f₁ is another arbitrary function. The final order we consider yields another forced problem:

O (η^{2}) : {\begin{cases} (Δ_{r, φ} + Ω_{0}^{2}) u_{2} = - \frac{2}{r^{2}} \frac{\partial^{2} u_{1}}{\partial φ \partial θ} - (\frac{1}{r^{2}} \frac{\partial^{2}}{\partial θ^{2}} + Ω_{2}^{2}) u_{0}, & (r, ϕ) \in S_{0}, \\ n \cdot \nabla_{r, φ} u_{2} = - n \cdot \hat{ϕ} \frac{1}{r} \frac{\partial u_{1}}{\partial θ}, & (r, ϕ) \in \cup \partial C_{j} \cap {\bar{S}}_{0} . \end{cases}

3.7

We derive a second compatibility condition at this order (appendix A(b)), this time yielding a second-order ODE satisfied by the modulation function f₀:

T \frac{d^{2} f_{0}}{d θ^{2}} + Ω_{2}^{2} f_{0} = 0,

3.8

where

T = \frac{\int_{S_{0}} (1 / r^{2}) (U_{0} (\partial U_{1} / \partial φ) - U_{1} (\partial U_{0} / \partial φ) + U_{0}^{2}) d S}{\int_{S_{0}} U_{0}^{2} d S + \int_{- π / N}^{π / N} (i / 2 Ω_{0}) r U_{0}^{2} |_{r = r_{+}} d φ} .

3.9

The asymptotic solutions constructed in this section will approximate solutions to the original system (2.1) provided that the value of r₊ has been chosen appropriately. A pair of linearly independent solutions to (3.8) is given by

f_{0} = \exp (\pm i Ω_{2} θ / \sqrt{T})

, and the cyclic continuity of the problem demands that

Ω_{2}^{2} / T = {(Δ m)}^{2} / η^{2}

for

Δ m \in Z

. This function governs the modulation of the periodic or antiperiodic standing wave solution; we have Δm=m if U₀ is periodic, and Δm=N/2−m if U₀ is antiperiodic. Using the ansatz (3.4), we derive an asymptotic expression for the complex frequencies Ω of the modulated quasi-modes, given the standing wave frequency Ω₀:

Ω = Ω_{0} + \frac{T {(Δ m)}^{2}}{2 Ω_{0}} + \dots

3.10

This expression is valid only for the smallest values of Δm=0,±1,±2,…, and ceases to be accurate when the scale-separation assumption no longer holds. For both of the systems we consider in this paper, only the most wide-angle modulated modes have Q-factors suggesting significant confinement, so the asymptotic method covers the regime of interest.

4. Geometric perturbation

The method outlined in §3 allows us to calculate asymptotic solutions to the rotationally symmetric system (2.1) that exhibit wide-angle Bloch modulation. In this section, we consider instead a case in which the geometry of the system is slightly perturbed such that this symmetry is broken. We could vary, for example, the length of or angle subtended by the spokes in the geometry of figure 2a, or the sizes or locations of the holes in figure 2b. We consider cases in which the size of the perturbation varies slowly from cell to cell; more precisely, the perturbation to the jth cell is proportional to a factor $g_{j} = \tilde{g} (2 π j / N)$ , where $\tilde{g} (ϕ)$ is a smooth, slowly varying function of the angular variable, and hence cyclic continuity demands that g is periodic such that g₀=g_N. If the geometric perturbation is sufficiently small, the result of the asymptotic procedure will be a Schrödinger-type ODE in place of (3.8), and the solutions of this will provide the underlying field distribution in the perturbed system. We consider a specific example to illustrate this.

(a) Angular perturbation in system of Neumann spokes

As a concrete example, let us consider introducing a variation in the angle subtended by each pair of adjacent spokes in the geometry of figure 2a. Suppose the angle between the jth and ( j+1)th spoke is given by

ϕ_{j + 1} - ϕ_{j} = \frac{2 π}{N} (1 + η^{2} g_{j}),

4.1

where

g_{j} = \tilde{g} (2 π j / N)

for a smooth, 2π-periodic function

\tilde{g} (ϕ)

, and

\sum_{j = 0}^{N - 1} g_{j} = 0

. Note that, unlike in the previous section, η≪1 must now be defined explicitly. In terms of the separated-scale variables,

\tilde{g} (ϕ)

is considered a function of the long-scale only, and we define

g (θ) = \tilde{g} (ϕ)

. We now introduce a scaled angular variable implicitly via ϕ=φ(1+η²g(θ)), where once again θ=ηϕ is the independent long-scale angular variable. In terms of the scaled variable φ, each cell subtends an angle of 2π/N, so we may consider an arbitrary cell as we did in the previous section. The azimuthal partial derivative operator is expanded in two scales using the chain rule as

\frac{\partial}{\partial ϕ} = {(1 + η^{2} g (θ))}^{- 1} \frac{\partial}{\partial φ} + η \frac{\partial}{\partial θ} = (1 - η^{2} g (θ)) \frac{\partial}{\partial φ} + η \frac{\partial}{\partial θ} + O (η^{4}) .

4.2

Having made this substitution, we seek solutions using the ansatz (3.4) as before, and find that the leading and first-order systems (3.5) and (3.6) are unchanged. The solutions U₀ and U₁ are carried forward to the second-order system, which is now given by

O (η^{2}) : {\begin{cases} (Δ_{r, φ} + Ω_{0}^{2}) u_{2} = - \frac{2}{r^{2}} \frac{\partial^{2} u_{1}}{\partial φ \partial θ} - (\frac{1}{r^{2}} \frac{\partial^{2}}{\partial θ^{2}} - \frac{2 g (θ)}{r^{2}} \frac{\partial^{2}}{\partial φ^{2}} + Ω_{2}^{2}) u_{0}, & (r, φ) \in S_{0}, \\ n \cdot \nabla_{r, φ} u_{2} = - n \cdot \hat{ϕ} \frac{1}{r} (\frac{\partial u_{1}}{\partial θ} - g (θ) \frac{\partial u_{0}}{\partial φ}), & (r, φ) \in \cup \partial C_{j} \cap {\bar{S}}_{0} . \end{cases}

4.3

The compatibility condition at this order is derived analogously to the one in (b), but now contains an extra term, resulting in an angular Schrödinger equation

T \frac{d^{2} f_{0}}{d θ^{2}} + (Ω_{2}^{2} - β g (θ)) f_{0} = 0,

4.4

where T is given by (3.9) and

β = \frac{\int_{S_{0}} (2 / r^{2}) U_{0} (d^{2} U_{0} / d φ^{2}) d S}{\int_{S_{0}} U_{0}^{2} d S + \int_{- π / N}^{π / N} (i / 2 Ω_{0}) r U_{0}^{2} |_{r = r_{+}} d φ} .

4.5

To illustrate the accuracy of the method developed in this section, in figure 9, we show an example where a system of N=30 spokes is perturbed as in (4.1), with the perturbation function given by $\tilde{g} (ϕ) = \cos (ϕ)$ , with η=2/N. In this case, the standing wave frequency is given by Ω₀=19.5237, and (4.4) yields a complex Mathieu equation, given explicitly by

(- 0.1474 - 0.1563 i) \frac{d^{2} f_{0}}{d θ^{2}} + (Ω_{2}^{2} + 647.36 \cos (15 θ)) f_{0} = 0.

4.6

Solving this using a spectral method [32], we find that the two most confined eigenvalues (those with the smallest imaginary parts) occur at

Ω_{2}^{2} = 534.61 - 47.28 i

and

Ω_{2}^{2} = 313.38 - 136.97 i

. These values are substituted into the ansatz (3.4) along with the standing wave frequency to give the predicted frequencies of the modes. Table 1 shows that these calculated values are very close to those calculated using the full FEM simulation.

Figure 9. The two most confined quasi-modes of a system of 30 spokes with an angular perturbation as described in (4.1), with $\tilde{g} (ϕ) = \cos (ϕ)$ with η=2/N. The real parts of the modulation functions f₀, calculated from equation (4.6), are shown in (a,b), with the real parts of the fields u, calculated by full FEM simulation, shown in (c,d).

Table 1.Frequency comparison of the two quasi-modes shown in figure 9, calculated first using the asymptotic theory and then full FEM simulation.
Collapse
mode from figure 9	Ω_Asymptotic	Ω_FEM
(c)	19.5845−0.0054i	19.5853−0.0053i
(d)	19.5594−0.0156i	19.5602−0.0155i

We note that the geometric perturbation above results in a significant reduction of the Q-factor compared with that of the standing wave solution in the unperturbed geometry. A similar phenomenon was observed in [7], and the asymptotic method presented here sheds light on this in the case of large numbers of inclusions: presuming that the Q-factor of the unperturbed standing wave is (formally) of order greater that N², the order 1/N² relative shift in real and imaginary parts of the frequency due to Ω₂ results in the Q-factor being reduced to $O (N^{2})$ . A detailed study of Q-factors in cyclic geometries will be presented in later work.

5. Concluding remarks

We have studied the eigenvalue problem for open systems with high degrees of rotational symmetry, governed by the planar Helmholtz equation. The WBM solutions have similarities with Rayleigh–Bloch surface waves, but fundamental difference can be identified: firstly, WBMs are quasi-modes, meaning they radiate energy to infinity and are hence characterized by complex frequencies. Secondly, the dispersion of the WBMs is discrete rather than continuous due to cyclic continuity. Thirdly, the standing wave WBMs can be either periodic or antiperiodic across one repeating cell, whereas the equivalent solutions for Rayleigh–Bloch surface waves are necessarily antiperiodic.

At frequencies close to those of the standing wave WBMs, solutions exist with wide-angle modulation, leading to dipolar and higher order radiation profiles. Analogous solutions are expected in other types of wave system, for example, those supporting in-plane elastic waves, bending waves in elastic plates and acoustic or electromagnetic waves in three dimensions. In the three-dimensional case, we expect WBMs to exist in finite structures with periodicity in the azimuthal direction. If similarly high degrees of confinement can be found for standing waves in these systems, the corresponding structures may be employed as alternatives to photonic or phononic crystal cavities in situations where confinement is desirable. Alternatively, a number of ring-like structures could be arranged on a Bravais lattice, in which case the WBM resonances could form the bases for new types of metamaterial or photonic crystal. Small perturbations of the sort considered in §4 may then be used to introduce directionality and other effects within these.

Author' contributions

Both of the authors have provided substantial contributions to the conception and design of the model, interpretation of the results and writing the article. The majority of the analysis was done by B.M. Both authors have given their final approval of the version to be published.

Competing interests

The authors of the paper have no competing interests.

Funding

This work was supported by EPSRC (UK) through Programme grant no. EP/L024926/1.

Acknowledgements

We are grateful to Alexander Movchan and Stewart Haslinger for useful conversations about multipole techniques and to the anonymous referees for comments that led to several improvements to the article.

Appendix A. Compatibility conditions in asymptotic procedure

(a) $O (η)$ compatibility condition

To derive the first compatibility condition, we multiply the first line of (3.5) by the unknown field u₁, and subtract the first line of (3.6) multiplied by u₀. The resulting equation is integrated over the cell $S_{0}$ , holding the long-scale variable θ constant:

\int_{S_{0}} (u_{1} Δ_{r, φ} u_{0} - u_{0} Δ_{r, φ} u_{1}) d S = \int_{S_{0}} (\frac{2}{r^{2}} u_{0} \frac{\partial^{2} u_{0}}{\partial φ \partial θ} + Ω_{1}^{2} u_{0}^{2}) d S .

A 1

With the help of a Green’s identity, the left-hand side can be written as a sum of line integrals. The contributions along the opposing edges of the cell cancel by (anti-)periodicity of the u_i’s, and after applying the boundary conditions in (3.5) and (3.6) on the inclusions, we are left with

\begin{aligned} {\int_{- π / N}^{π / N} (u_{1} \frac{\partial u_{0}}{\partial r} - u_{0} \frac{\partial u_{1}}{\partial r}) r |}_{r = r_{+}} d φ + \int_{\cup \partial C_{j} \cap {\bar{S}}_{0}} n \cdot \hat{ϕ} \frac{1}{r} u_{0} \frac{\partial u_{0}}{\partial θ} d l \\ = \int_{S_{0}} (\frac{2}{r^{2}} u_{0} \frac{\partial^{2} u_{0}}{\partial φ \partial θ} + Ω_{1}^{2} u_{0}^{2}) d S . \end{aligned}

A 2

We substitute in the separated-scale leading order solution u₀(r,φ,θ)=f₀(θ)U₀(r,φ) and rearrange, which yields

\begin{aligned} {\int_{- π / N}^{π / N} (u_{1} \frac{\partial U_{0}}{\partial r} - U_{0} \frac{\partial u_{1}}{\partial r}) r |}_{r = r_{+}} d φ - Ω_{1}^{2} f_{0} \int_{S_{0}} U_{0}^{2} d S \\ = \frac{d f_{0}}{d θ} (\int_{S_{0}} \frac{2}{r^{2}} U_{0} \frac{\partial U_{0}}{\partial φ} d S - \int_{\cup \partial C_{j} \cap {\bar{S}}_{0}} n \cdot \hat{ϕ} \frac{1}{r} U_{0}^{2} d l) . \end{aligned}

A 3

Integrating by parts, we find that the terms on the right-hand side sum to zero, so we are left with

{\int_{- π / N}^{π / N} (u_{1} \frac{\partial U_{0}}{\partial r} - U_{0} \frac{\partial u_{1}}{\partial r}) r |}_{r = r_{+}} d φ - Ω_{1}^{2} f_{0} \int_{S_{0}} U_{0}^{2} d S = 0.

A 4

At r=r₊, the field u satisfies the outgoing wave condition ∂u/∂r=(iΩ−1/2r)u. Expanding this to

O (η)

using the ansatz (3.4) yields

\frac{\partial u}{\partial r} = (i Ω_{0} u_{0} - \frac{u_{0}}{2 r}) + η (i Ω_{0} u_{1} + i \frac{Ω_{1}^{2}}{2 Ω_{0}} u_{0} - \frac{u_{1}}{2 r}) + O (η^{2}),

A 5

where the terms in brackets are identified as ∂u_j/∂r for j=0,1 respectively. Substituting these expressions for the derivatives into (A 4) yields the result

C Ω_{1}^{2} \equiv (\int_{S_{0}} U_{0}^{2} d S + {\int_{- π / N}^{π / N} \frac{i}{2 Ω_{0}} r U_{0}^{2} |}_{r = r_{+}} d φ) Ω_{1}^{2} = 0.

A 6

For non-trivial solutions, the quantity C is non-zero, so we deduce that Ω₁=0.

(b) $O (η^{2})$ compatibility condition

We proceed by the same method as in appendix A(a), multiplying the first line of (3.5) by the unknown field u₂, and subtracting the first line of (3.7) multiplied by u₀. The resulting equation is integrated over the cell $S_{0}$ . We again use a Green’s identity on the left-hand side, and after applying the boundary conditions and substituting the separated-scale expressions for u₀ and u₁ are left with the equation

\begin{aligned} {\int_{- π / N}^{π / N} (u_{2} \frac{\partial U_{0}}{\partial r} - U_{0} \frac{\partial u_{2}}{\partial r}) r |}_{r = r_{+}} d φ - Ω_{2}^{2} f_{0} \int_{S_{0}} U_{0}^{2} d S \\ = \frac{d^{2} f_{0}}{d θ^{2}} \int_{S_{0}} \frac{1}{r^{2}} (U_{0} \frac{\partial U_{1}}{\partial φ} - U_{1} \frac{\partial U_{0}}{\partial ϕ} + U_{0}^{2}) d S . \end{aligned}

A 7

Expanding the outgoing wave condition again, this time to

O (η^{2})

, yields

\frac{\partial u}{\partial r} = (i Ω_{0} u_{0} - \frac{u_{0}}{2 r}) + η (i Ω_{0} u_{1} - \frac{u_{1}}{2 r}) + η^{2} (i Ω_{0} u_{2} + i \frac{Ω_{2}^{2}}{2 Ω_{0}} u_{0} - \frac{u_{2}}{2 r}) + O (η^{3})

A 8

and substituting these expressions for the derivatives into (A 7) yields the result

\begin{aligned} C Ω_{2}^{2} f_{0} & \equiv ({\int_{S_{0}} U_{0}^{2} d S + \int_{- π / N}^{π / N} \frac{i}{2 Ω_{0}} r U_{0}^{2} |}_{r = r_{+}} d φ) Ω_{2}^{2} f_{0} \\ = - \frac{d^{2} f_{0}}{d θ^{2}} \int_{S_{0}} \frac{1}{r^{2}} (U_{0} \frac{\partial U_{1}}{\partial φ} - U_{1} \frac{\partial U_{0}}{\partial φ} + U_{0}^{2}) d S, \end{aligned}

A 9

which is rearranged to give (3.8). We stress that this equation governs the angular modulation of solutions to the truncated-cell problem. For a sufficiently large value of r₊, these approximate solutions to the full planar problem (2.1), which are the ones of interest. We note, however, that if r₊ is chosen to be too large, the asymptotic hierarchy will no longer hold, causing the method to fail. To understand why this is the case, we note that the converged solutions satisfy

U_{j} \propto e^{i Ω_{0} r} / \sqrt{r}

for j=0,1 as r→r₊. Since Ω₀ has a negative imaginary part, the outgoing waves grow exponentially as

r_{+} \to \infty

, causing the hierarchy to break down and the integrals in (A 9) to diverge. Since the amplitudes of these waves are minimum at r=1/|2 Im (Ω₀)|, the hierarchy will remain intact if we choose r₊<1/|2 Im (Ω₀)| and this issue will be avoided.

Assuming |Im (Ω₀)|≪1, there is freedom to choose r₊ in the regime 1≪r₊<1/|2Im (Ω₀)|. It will be useful to have a bound on the error associated with making a particular choice, compared with the ‘best’ case r₊=1/|2Im (Ω₀)|, which may be excessively large and inconvenient for numerical methods such as FEM. Assuming the solution is sufficiently converged, the quantity $C \in C$ appearing on the left-hand side of (A 9) is independent of r₊, as shown by considering the derivative:

\frac{d C}{d r_{+}} = {\int_{- π / N}^{π / N} (r U_{0}^{2} + \frac{i}{2 Ω_{0}} \frac{d}{d r} (r U_{0}^{2})) |}_{r = r_{+}} d φ .

A 10

Since ∂U₀/∂r=(iΩ₀−1/2r)U₀ at r=r₊, we have

{\frac{d}{d r} (r U_{0}^{2}) |}_{r = r_{+}} = {(U_{0}^{2} + 2 r U_{0}^{2} (i Ω_{0} - \frac{1}{2 r})) |}_{r = r_{+}} = 2 i Ω_{0} r U_{0}^{2} |_{r = r_{+}},

A 11

which is substituted into (A 10) to yield the result dC/dr₊=0, and hence the left-hand side of (A 9) is independent of r₊.

Next, we consider the right-hand side of (A 9). For converged solutions, the magnitude of the extra contribution to the integral in moving the boundary from r₊ to r₊′ is estimated by

| A \int_{r_{+}}^{r_{+}^{'}} \frac{1}{r^{2}} {(\frac{e^{i Ω_{0} r}}{\sqrt{r}})}^{2} r d r | < max_{r_{+} < q < r_{+}^{'}} | A \frac{1}{r_{+}^{2}} \int_{r_{+}}^{q} e^{2 i Ω_{0} r} d r | < | \frac{A}{Ω_{0} r_{+}^{2}} | e^{2 | Im (Ω_{0}) | r_{+}^{'}},

A 12

where

A = \int_{- π / N}^{π / N} e^{- 2 i Ω_{0} r} (U_{0} (\partial U_{1} / \partial ϕ) - U_{1} (\partial U_{0} / \partial ϕ) + U_{0}^{2}) r |_{r = r_{+}} d ϕ

Since we have also demanded that r₊′<1/|2Im (Ω₀)|, the variation in T is bounded by |Ae/(CΩ₀r²₊)|. A value of r₊ should be chosen such that this value is significantly smaller than the higher-order terms neglected in the asymptotic expansion.

Footnotes

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

References

1
Rayleigh L. 1910The problem of the whispering gallery. Philos. Mag. 20, 1001–1004. (doi:10.1080/14786441008636993) Crossref, Google Scholar
2
Matsko AB, Savchenkov AA, Strekalov D, Ilchenko VS, Maleki L. 2005Review of applications of whispering-gallery mode resonators in photonics and nonlinear optics. IPN Progress Report. Google Scholar
3
McCall SL, Levi AFJ, Slusher RE, Pearton SJ, Logan RA. 1992Whispering gallery mode microdisk lasers. App. Phys. Lett. 60, 289. (doi:10.1063/1.106688) Crossref, ISI, Google Scholar
4
Vollmer F, Arnold S. 2008Whispering-gallery-mode biosensing: label-free detection down to single molecules. Nat. Methods 5, 591–596. (doi:10.1038/nmeth.1221) Crossref, PubMed, ISI, Google Scholar
5
Ilchenko VS, Savchenkov AA, Matsko AB, Maleki L. 2004Nonlinear optics and crystalline whispering gallery mode cavities. Phys. Rev. Lett. 92, 043903. (doi:10.1117/12.538445) Crossref, PubMed, ISI, Google Scholar
6
McIver P, Linton CM, McIver M. 1998Construction of trapped modes for wave guides and diffraction gratings. Proc. R. Soc. Lond. A 454, 2593–2616. (doi:10.1098/rspa.1998.0272) Link, ISI, Google Scholar
7
Evans D, Porter R. 1997Near trapping of waves by circular arrays of vertical cylinders. Appl. Ocean Res. 19, 83–99. (doi:10.1016/S0141-1187(97)00015-1) Crossref, ISI, Google Scholar
8
Chapman SJ, Hewett DP, Trefethen LN. 2015Mathematics of the Faraday cage. SIAM Rev. 57, 398–417. (doi:10.1137/140984452) Crossref, ISI, Google Scholar
9
Martin PA. 2014On acoustic and electric Faraday cages. Proc. R. Soc. A 470, 20140344. (doi:10.1098/rspa.2014.0344) Link, Google Scholar
10
Antonakakis T, Craster RV, Guenneau S, Skelton EA. 2014An asymptotic theory for waves guided by diffraction gratings or along microstructured surfaces. Proc. R. Soc. A 470, 20130467. (doi:10.1098/rspa.2013.0467) Link, Google Scholar
11
Evans DV, Linton CM. 1993Edge waves along periodic coastlines. Q. J. Appl. Math. 46, 643–656. (doi:10.1093/qjmam/46.4.643) Crossref, ISI, Google Scholar
12
Porter R, Evans DV. 2005Embedded Rayleigh-Bloch surface waves along periodic rectangular arrays. Wave Motion 43, 29–50. (doi:10.1016/j.wavemoti.2005.05.005) Crossref, ISI, Google Scholar
13
Porter R, Evans DV. 1999Rayleigh-Bloch surface waves along periodic gratings and their connection with trapped modes in waveguides. J. Fluid Mech. 386, 233–258. (doi:10.1017/S0022112099004425) Crossref, ISI, Google Scholar
14
Linton CM, McIver M. 2002The existence of Rayleigh-Bloch surface waves. J. Fluid Mech. 470, 85–90. (doi:10.1017/S0022112002002227) Crossref, ISI, Google Scholar
15
Bonnet-Bendhia AS, Starling F. 1994Guided waves by electromagnetic gratings and non- uniqueness examples for the diffraction problem. Math. Meth. Appl. Sci. 17, 305–338. (doi:10.1002/mma.1670170502) Crossref, ISI, Google Scholar
16
Linton CM, Zalipaev V, Thompson I. 2013Electromagnetic guided waves on linear arrays of spheres. Wave Motion 50, 29–40. (doi:10.1016/j.wavemoti.2012.06.002) Crossref, ISI, Google Scholar
17
Fernandez-Dominguez AI, Garcia-Vidal F, Martin-Moreno L. 2011Surface electromagnetic waves on structured perfectly conducting surfaces. In Structured surfaces as optical metamaterials (ed. AA Maradudin). Cambridge, UK: Cambridge University Press. Google Scholar
18
Pendry JB, Martin-Moreno L, Garcia-Vidal FJ. 2004Mimicking surface plasmons with structured surfaces. Science 305, 847–848. (doi:10.1126/science.1098999) Crossref, PubMed, ISI, Google Scholar
19
Huidubro PA, Xiaopeng S, Cuerda J, Moreno E, Martin-Moreno L, Garcia-Vidal FJ, Cui TJ, Pendry JB. 2014Magnetic localized surface plasmons. Phys. Rev. X 4, 021003. (doi:10.1103/PhysRevX.4.021003) ISI, Google Scholar
20
Hurd RA. 1954The propagation of an electromagnetic wave along an infinite corrugated surface. Can. J. Phys. 32, 727–734. (doi:10.1139/p54-079) Crossref, ISI, Google Scholar
21
Sengupta D. 1959On the phase velocity of wave propagation along an infinite yagi structure. IRE Trans. Antennas Propag. 7, 234–239. (doi:10.1109/TAP.1959.1144682) Crossref, Google Scholar
22
Makwana M, Antonakakis A, Maling B, Guenneau S, Craster RV. 2016Wave mechanics in media pinned at Bravais lattice points. SIAM J. Appl. Math. 76, 1–26. (doi:10.1137/15M1020976) Crossref, ISI, Google Scholar
23
Craster RV, Kaplunov J, Pichugin AV. 2010High-frequency homogenization for periodic media. Proc. R. Soc. A 466, 2341–2362. (doi:10.1098/rspa.2009.0612) Link, Google Scholar
24
Brillouin L. 1946Wave propagation in periodic structures. New York, NY: Dover. Google Scholar
25
Olver FWJ. 1952Some new asymptotic expansions for Bessel functions of large orders. Proc. Camb. Philos. 48, 414–427. (doi:10.1017/S030500410002781X) Crossref, Google Scholar
26
Movchan AB, Movchan NV, Poulton CG. 2002Asymptotic models of fields in dilute and densely packed composites. London, UK: Imperial College Press. Crossref, Google Scholar
27
Linton CM, McIver P. 2001Handbook of mathematical techniques for wave/structure interactions. London, UK: Chapman & Hall CRC. Crossref, Google Scholar
28
Wijngaard W. 1973Guided normal modes of two parallel circular dielectric rods. J. Opt. Soc. Am. 63, 944–949. (doi:10.1364/JOSA.63.000944) Crossref, ISI, Google Scholar
29
Abramowitz M, Stegun IA. 1964Handbook of mathemtical functions. New York, NY: Dover. Google Scholar
30
Colquitt DJ, Craster RV, Antonakakis T, Guenneau S. 2014Rayleigh–Bloch waves along elastic diffraction gratings. Proc. R. Soc. A 471, 20140465. (doi:10.1098/rspa.2014.0465) Link, Google Scholar
31
COMSOL Ltd. 2014Comsol multiphysics 5.0. Google Scholar
32
Trefethen LN. 2000Spectral methods in MATLAB. Philadelphia, PA: SIAM. Crossref, Google Scholar