## Abstract

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 *ω*)<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 1*a*,*b*. In the case with spokes, however, there is an additional mode, shown in figure 1*c*, that is created by the periodicity of the geometry; this second mode persists even if the materials are identical, as shown in figure 1*d*. 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.

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 *N*. Note that *N*, and a simple case to consider would be that of **n**, and *r*=|**x**|. In particular, we consider the two geometries shown in figure 2 so that concrete solutions can be found.

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 *u*(*r*,*ϕ*), and since the eigenvalue is simple, the two must be identical up to multiplication by a complex number, i.e. *N* times, the cyclic continuity of the problem requires that *α*=e^{2mπi/N} is one of the *N*th roots of unity. The rotational analogue of Bloch’s theorem is thus
*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 2*a*.

### (a) Neumann spokes: matched eigenfunction expansion method

We consider the infinite wedge-shaped repeating cell as in figure 2*a*. 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 *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 *ϕ*. Using the orthogonality of the cosine functions, along with the relation

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 *M*+1 equations with −*M*≤*n*≤*M*, 0≤*p*≤2*M*, 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 *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*)=(−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 *m*, there will be a spectrum **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 2*b*, 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 *j*th term in the first summation is given in terms of coordinates centred at the *j*th hole, and hence *u* is expressed as a function of 2*N* (dependent) variables. This expansion is valid everywhere in *r*_{0j} is the distance between the centres of the 0th and *j*th 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, *M*+1 equations with −*M*≤*μ*,*ν*≤*M*. The corresponding matrix equation **PA**=**0** has non-trivial solutions iff *m*, there will be a spectrum **A**.

## 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 *r*_{+}≫1. We impose the outgoing wave condition ∂*u*/∂*r*=(*iΩ*−1/2*r*)*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:
*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
*et al.* [23], we seek solutions in the form of an asymptotic ansatz:
*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

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*_{0}(*r*,*φ*), multiplied by a coefficient that generally depends on the independent long-scale angular coordinate *θ*, so we have *u*_{0}(*r*,*φ*,*θ*)=*f*_{0}(*θ*)*U*_{0}(*r*,*φ*) for an arbitrary function *f*_{0}. 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 *Ω*_{0} with those calculated using the expansion methods in §2. Since we are interested in spatially confined quasi-modes, it is assumed that *Ω*_{0} has only a small imaginary part. Specifically, the condition that *r*_{+}<1/|2Im (*Ω*_{0})| 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:
*Ω*_{1}=0 (appendix A(a)). With this in place, the forced system is solved, yielding a solution of the form *u*_{1}(*r*,*φ*,*θ*)=*f*_{1}(*θ*)*U*_{0}(*r*,*φ*)+*f*_{0}′(*θ*)*U*_{1}(*r*,*φ*), where *f*_{1} is another arbitrary function. The final order we consider yields another forced problem:
*f*_{0}:
*r*_{+} has been chosen appropriately. A pair of linearly independent solutions to (3.8) is given by *m*=*m* if *U*_{0} is periodic, and Δ*m*=*N*/2−*m* if *U*_{0} 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}:
*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 2*a*, or the sizes or locations of the holes in figure 2*b*. We consider cases in which the size of the perturbation varies slowly from cell to cell; more precisely, the perturbation to the *j*th cell is proportional to a factor *g* is periodic such that *g*_{0}=*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 2*a*. Suppose the angle between the *j*th and ( *j*+1)th spoke is given by
*π*-periodic function *η*≪1 must now be defined explicitly. In terms of the separated-scale variables, *ϕ*=*φ*(1+*η*^{2}*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
*U*_{0} and *U*_{1} are carried forward to the second-order system, which is now given by
*T* is given by (3.9) and

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 *η*=2/*N*. In this case, the standing wave frequency is given by *Ω*_{0}=19.5237, and (4.4) yields a complex Mathieu equation, given explicitly by

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*^{2}, the order 1/*N*^{2} relative shift in real and imaginary parts of the frequency due to *Ω*_{2} results in the *Q*-factor being reduced to *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) **

To derive the first compatibility condition, we multiply the first line of (3.5) by the unknown field *u*_{1}, and subtract the first line of (3.6) multiplied by *u*_{0}. The resulting equation is integrated over the cell *θ* constant:
*u*_{i}’s, and after applying the boundary conditions in (3.5) and (3.6) on the inclusions, we are left with
*u*_{0}(*r*,*φ*,*θ*)=*f*_{0}(*θ*)*U*_{0}(*r*,*φ*) and rearrange, which yields
*r*=*r*_{+}, the field *u* satisfies the outgoing wave condition ∂*u*/∂*r*=(*iΩ*−1/2*r*)*u*. Expanding this to *u*_{j}/∂*r* for *j*=0,1 respectively. Substituting these expressions for the derivatives into (A 4) yields the result
*C* is non-zero, so we deduce that *Ω*_{1}=0.

**(b) **

We proceed by the same method as in appendix A(a), multiplying the first line of (3.5) by the unknown field *u*_{2}, and subtracting the first line of (3.7) multiplied by *u*_{0}. The resulting equation is integrated over the cell *u*_{0} and *u*_{1} are left with the equation
*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 *j*=0,1 as *r*→*r*_{+}. Since *Ω*_{0} has a negative imaginary part, the outgoing waves grow exponentially as *r*=1/|2 Im (*Ω*_{0})|, the hierarchy will remain intact if we choose *r*_{+}<1/|2 Im (*Ω*_{0})| and this issue will be avoided.

Assuming |Im (*Ω*_{0})|≪1, there is freedom to choose *r*_{+} in the regime 1≪*r*_{+}<1/|2Im (*Ω*_{0})|. 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 (*Ω*_{0})|, which may be excessively large and inconvenient for numerical methods such as FEM. Assuming the solution is sufficiently converged, the quantity *r*_{+}, as shown by considering the derivative:
*U*_{0}/∂*r*=(*iΩ*_{0}−1/2*r*)*U*_{0} at *r*=*r*_{+}, we have
*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

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

- Received February 11, 2016.
- Accepted June 17, 2016.

- © 2016 The Authors.

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.