Whispering Bloch modes

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.


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  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.
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.
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 wideangle 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.

Prototype system
We pose an eigenvalue problem for the planar Helmholtz equation with Neumann inclusions, seeking u : Here C j for j = 1, . . . ,Ñ refers to a set of bounded inclusions arranged to have rotational symmetry of order N. Note thatÑ may be greater or less than N, and a simple case to consider would be that ofÑ = 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. 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 operatorR such that for any suitably defined function,Rf (r, φ) = f (r, φ + 2π/N). The rotational symmetry of (2.1) means that the fieldRu(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.Ru(r, φ) = αu(r, φ) for α ∈ C. If we perform the same transformation N times, the cyclic continuity of the problem requires that α N u(r, φ) =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 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 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(pNφ/2) for p ∈ Z and integrate over the azimuthal variable φ. Using the orthogonality of the cosine functions, (b) (a) 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.
along with the relation we are able to eliminate the coefficientsÂ n , leaving an infinite system of equations: and μ (2) 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  B n ∼ O(n −3/2 ) as n → ∞. 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 pn =P pn /(J pN/2 (Ω)H m+nN (Ω)), which satisfy 0 = n P pn B n ∀ p ∈ Z for the same values of Ω ∈ 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]: (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  Figure 5. Geometry of the cell S j . In order for the holes C j and C j±1 to be non-overlapping, we require that a < sin(π/N). 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

Re
where A 0 μ J μ (Ωa) + E 0 μ H μ (Ωa) = 0 ∀ μ ∈ 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 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 \ ∪ C i , and automatically satisfies the outgoing radiation condition as r → ∞.
Quasi-periodicity implies that E j ν = E 0 ν e 2jmπi/N ∀ ν ∈ Z. We use Graf's addition theorem [29] to express every term in (2.10) in the coordinate system of the 0th hole: 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,

Asymptotic analysis
The existence of standing wave WBMs, which are periodic or antiperiodic across a single wedgeshaped 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 wedgeshaped 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 ∈ (0, r + ], φ ∈ [−π/N, π/N]}\ ∪ 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: 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 and Following the approach of Craster et al. [23], we seek solutions in the form of an asymptotic ansatz: u(r, ϕ, θ) = u 0 (r, ϕ, θ) + ηu 1 (r, ϕ, θ) + η 2 u 2 (r, ϕ, θ) + · · · and 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 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: At this point, we derive a compatibility condition that implies Ω 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: 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 0 : T 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. of the periodic or antiperiodic standing wave solution; we have 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 : 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.

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 =g(2π j/N), whereg(φ) is a smooth, slowly varying function of the angular variable, and hence cyclic continuity demands that 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 2a. Suppose the angle between the jth and ( j + 1)th spoke is given by where g j =g(2π j/N) for a smooth, 2π -periodic functiong(φ), and N−1 j=0 g j = 0. Note that, unlike in the previous section, η 1 must now be defined explicitly. In terms of the separated-scale variables,g(φ) is considered a function of the long-scale only, and we define g(θ) =g(φ). We now introduce a scaled angular variable implicitly via φ = ϕ(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 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 0 and U 1 are carried forward to the second-order system, which is now given by 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 where 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 byg(φ) = cos(φ), with η = 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 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.28i and Ω 2 2 = 313.38 − 136.97i. 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.
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 O(N 2 ). A detailed study of Q-factors in cyclic geometries will be presented in later work.

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.
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 ∝ e iΩ 0 r / √ r for j = 0, 1 as r → r + . Since Ω 0 has a negative imaginary part, the outgoing waves grow exponentially as r + → ∞, 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 (Ω 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 C ∈ C appearing on the left-hand side of (A 9) is independent of r + , as shown by considering the derivative: dC dr + = 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 where A = π/N −π/N e −2iΩ 0 r (U 0 (∂U 1 /∂φ) − U 1 (∂U 0 /∂φ) + U 2 0 )r| r=r + dφ. 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.