#  Two-dimensional Helmholtz resonator arrays. Part II. Matched asymptotic
expansions for specially-scaled resonators

M. J. A. Smith  Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Wilberforce Road, CB3 0WA, UK I. D. Abrahams
Department of Applied Mathematics and Theoretical Physics, University of
Cambridge, Wilberforce Road, CB3 0WA, UK

###### Abstract

We present a solution method which combines the method of matched asymptotics
with the method of multipole expansions to determine the band structure of
cylindrical Helmholtz resonators arrays in two dimensions. The resonator
geometry is considered in the limit as the wall thickness becomes very large
compared with the aperture width (the specially-scaled limit). In this regime,
the existing treatment in Part I, with updated parameters, is found to return
spurious spectral behaviour. We derive a regularised system which overcomes
this issue and also derive compact asymptotic descriptions for the low-
frequency dispersion equation in this setting. In the specially-scaled limit,
our asymptotic dispersion equation not only recovers the first band surface
but also extends to high, but still subwavelength, frequencies. A
homogenisation treatment is outlined for describing the effective bulk modulus
and effective density tensor of the resonator array for all wall thicknesses.
We demonstrate that specially-scaled resonators are able to achieve
exceptionally low Helmholtz resonant frequencies, and present closed-form
expressions for determining these explicitly. We anticipate that the
analytical expressions and the formulation outlined here may prove useful in
industrial and other applications.

Submitted Manuscript

##  1 Introduction

In Part I of this study, we outlined a matched asymptotic-multipole treatment
for determining the band structure of thin- and moderately thick-walled
Helmholtz resonator arrays. However, this formulation implicitly assumed that
the wall thickness was not too large as compared with the aperture width; an
assumption that prevents us from achieving very low Helmholtz resonance
frequencies, or equivalently, very low first-band gaps. Here we consider an
important extension to the results derived in Part I by examining arrays of
specially-scaled thick-walled resonators (with aspect ratios chosen to achieve
a low frequency resonance – hence the use of the terminology specially-scaled
resonators), see Figure 1. We also discuss a homogenisation procedure for all
wall thickness configurations.

On the topic of homogenisation, we remark that the literature on two- and
three-dimensional Helmholtz resonator arrays is extensive, particularly since
the resonator can exhibit (1) singular behaviour with respect to its geometry,
as well as (2) singular behaviour with respect to a constitutive quantity
(e.g., high-contrast conductivity). We do not attempt to conduct an exhaustive
literature review here but will instead highlight key works of interest. To
clarify by means of example, a geometrically singular medium could include an
array where the aperture width of each resonator contracts much more rapidly
than some other geometric parameter, as wavelengths become long relative to
the period of the unit cell (see works in acoustics [1, 2] and in
electromagnetism [3]). Interestingly, existing work in acoustics [1, 2]
suggest that for the specially-scaled Helmholtz resonator array problem, the
effective density is not frequency dependent whereas the effective bulk
modulus is indeed frequency dependent and takes negative values for a fixed
frequency interval. This very same behaviour is observed in two-scale
asymptotic treatments of high contrast arrays of cylinders for the Helmholtz
equation [4, 5]. The observation that frequency dependence emerges only in the
effective bulk modulus is consistent with the assertion in [6] that a dynamic
compressibility is observed near the Helmholtz resonance frequency for
Helmholtz resonator arrays, and not a dynamic density response.

For those interested in the history of the Helmholtz resonator design, it
would seem that the systematic investigation of cavity resonators emerged in
the mid-19th century, after tone generation was observed in heated glass
flasks with long necks [7]. Soon thereafter, mathematical explanations were
given by Helmholtz [8] and Lord Rayleigh [9], the first of these giving rise
to the nomenclature of Helmholtz resonators that is used to the present day.
The study of two-dimensional resonators seemingly have their origins in the
study of loop antennas which emerged almost a century later [10], from which
the nomenclature of split-ring resonator likewise emerged for the design
considered here, as understood by many in the community. For reference,
alternative nomenclatures include loop-gap resonator or split-tube resonator
within the electrical engineering literature.

The outline of this paper is as follows. First, we briefly restate the
governing equations offered in Part I in the following section. We then
determine the leading-order outer and inner solutions within the neck region
in Section 2, where asymptotic matching is also conducted. In Section 3 we
construct the regularised system for specially-scaled resonator arrays before
constructing asymptotic dispersion equations in Section 4. In Section 5 we
present closed-form representations for the Helmholtz resonance/cutoff
frequency. Numerical examples are then considered in Section 7, and finally an
extended discussion is given in Section 8, which highlights the differences of
the model results presented here to those presented in Part I.

###  1.1 Governing wave equation

We consider wave propagation in an acoustic medium satisfying the scalar wave
equation

$(\partial_{\overline{x}}^{2}+\partial_{\overline{y}}^{2})\,\overline{\phi}+k^{2}\overline{\phi}=0,$
(1)

where $\overline{x}$ denotes dimensional Cartesian coordinates,
$k^{2}=\rho\omega^{2}/B$ is the square of the wave number, $B$ is the Bulk
modulus, $\rho$ is the mass density of the medium, and $\bar{\phi}$ is the
velocity potential (see Part I for further details). Within this acoustic
medium we immerse a two-dimensional square array of rigid resonators as shown
in Figure 1, spaced a distance $\overline{d}$ apart, satisfying Neumann
conditions on the wall edges, and Bloch conditions throughout the cell.
Implicitly we examine time-harmonic solutions of the form
$\mathrm{exp}(-\mathrm{i}\omega t)$ where $\omega$ is the angular frequency,
but this factor is suppressed throughout for ease of exposition.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

Figure 1:  LABEL:sub@fig:ucell1 Representative fundamental unit cell in non-
dimensional coordinates for a square array of Helmholtz resonators with large
aspect ratios (i.e., specially-scaled Helmholtz resonators), where $a$ and $b$
denote the inner and outer radii, respectively, and $d$ is the periodicity;
LABEL:sub@fig:ucell2 Close-up of neck region showing central aperture angle
$\theta_{0}$ and aperture half-width angle $\theta_{\mathrm{ap}}$.

##  2 Helmholtz resonators in the specially-scaled limit

We shall now employ the method of matched asymptotic expansions [11, 12] as
outlined in Part I, where through nondimensional rescaling we introduced inner
and outer regions of the unit cell. However, in contrast to the treatment of
thick-walled resonators in Part I, where we partitioned the unit cell for the
outer problem into two domains, we now partition the unit cell for the outer
problem into three domains: the interior, neck, and exterior regions, and
consider the solution in each domain. The outer solutions for the interior and
exterior regions are identical to those presented in Part I, whose results we
restate for reference below. For the outer domains we use coordinates scaled
on the wave number

$x=k\overline{x},\quad\mbox{and}\quad y=k\overline{y},$ (2)

with the inner and outer radii of the resonator, $\overline{a}$ and
$\overline{b}$, scaled as $a=k\overline{a}$ and $b=k\overline{b}$
respectively, the lattice period $\overline{d}$ is scaled as
$d=k\overline{d}$, the nondimensional cylinder thickness is
$2\mathbcal{m}=b-a$, and the aperture half-width $\overline{\ell}$ is scaled
as $\varepsilon=k\ell$. In all that follows we take the asymptotic limit
$\varepsilon\rightarrow 0$. We now examine the outer solution in the neck
region in detail.

###  2.1 Outer solution in the neck region

Note that from the unit cell configuration shown in Figure 1 we first rotate
and translate the lattice as
$(\widetilde{x},\widetilde{y})\mapsto(x\sin\theta_{0}-y\cos\theta_{0},x\cos\theta_{0}+y\sin\theta_{0}-b)$,
where $\theta_{0}$ is the central aperture angle, so that the exterior mouth
of the resonator is located at the origin in $(\widetilde{x},\widetilde{y})$
coordinates. In this new coordinate frame we solve Helmholtz’s equation

$(\partial_{\widetilde{x}}^{2}+\partial_{\widetilde{y}}^{2}+1)\phi_{\mathrm{neck}}=0,$
(3)

in the neck region of the resonator
$S_{N}=\left\\{(\widetilde{x},\widetilde{y}):(-\varepsilon,\varepsilon)\times(-2\mathbcal{m},0)\right\\}$
subject to the Neumann boundary conditions

$\frac{\partial\phi_{\mathrm{neck}}}{\partial\widetilde{x}}\bigg{|}_{\widetilde{x}=\pm\varepsilon}=0,$
(4)

as illustrated in Figure 2, in the limit as $\varepsilon\rightarrow 0$. Away
from the aperture mouths, (3) and (4) admit the general solution

$\phi_{\mathrm{neck}}=\sum_{n=0}^{\infty}\left\\{p_{n}\mathrm{e}^{\mathrm{i}\lambda_{n}\widetilde{y}}+q_{n}\mathrm{e}^{-\mathrm{i}\lambda_{n}\widetilde{y}}\right\\}\cos\left(\frac{n\pi(\widetilde{x}+\varepsilon)}{2\varepsilon}\right),\qquad$
(5)

where

$\lambda_{n}=\begin{cases}\sqrt{1-\left(\dfrac{n\pi}{2\varepsilon}\right)^{2}},&\dfrac{n\pi}{2\varepsilon}<1,\\\
\mathrm{i}\sqrt{\left(\dfrac{n\pi}{2\varepsilon}\right)^{2}-1},&\dfrac{n\pi}{2\varepsilon}>1.\end{cases}$
(6)

In the closing aperture limit, the dominant contribution comes from the $n=0$
term and so the solution takes the form

$\lim_{\varepsilon\rightarrow 0}\phi_{\mathrm{neck}}\sim
p_{0}\mathrm{e}^{\mathrm{i}\widetilde{y}}+q_{0}\mathrm{e}^{-\mathrm{i}\widetilde{y}}.$
(7)

Accordingly the outer solution asymptotics near the entrance and exit to the
neck are given by

$\displaystyle\lim_{\widetilde{x}\rightarrow 0}\lim_{\widetilde{y}\rightarrow
0}\phi_{\mathrm{neck}}$
$\displaystyle\sim(p_{0}+q_{0})+\mathrm{i}\widetilde{y}(p_{0}-q_{0}),$ (8a)
$\displaystyle\lim_{\check{x}\rightarrow 0}\lim_{\check{y}\rightarrow
0}\phi_{\mathrm{neck}}$
$\displaystyle\sim\left(p_{0}\mathrm{e}^{-2\mathrm{i}\mathbcal{m}}+q_{0}\mathrm{e}^{2\mathrm{i}\mathbcal{m}}\right)+\mathrm{i}\check{y}(p_{0}\mathrm{e}^{-2\mathrm{i}\mathbcal{m}}-q_{0}\mathrm{e}^{2\mathrm{i}\mathbcal{m}}),$
(8b)

respectively, where we express (8b) in terms of the shifted origin
$(\check{x},\check{y})=(\widetilde{x},\widetilde{y}+2\mathbcal{m})$. We note
that for the inner solutions that follow, we use the same rotated and
translated frame $(\widetilde{x},\widetilde{y})$, and so we do not need to
express the representations (8) above in terms of $(x,y)$.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

Figure 2:  LABEL:sub@fig:outerneckthickres1 Outer problem inside the thick-
walled resonator neck of width $2\varepsilon$ and length $2\mathbcal{m}$ (not
to scale) where $(\widetilde{x},\widetilde{y})$ denotes the centre of the
exterior mouth and $(\check{x},\check{y})$ the centre of the interior mouth,
LABEL:sub@fig:outerneckthickres2(i) Inner problem geometry at the exterior
neck entrance, and LABEL:sub@fig:outerneckthickres2(ii) Inner problem geometry
after applying the Schwartz–Christoffel mapping (11); the capital letters
$A,\ldots,D$ and $A^{\prime},\ldots,D^{\prime}$ denote points of
correspondence between the two complex planes

###  2.2 Outer solution in the interior and exterior regions

Next, we restate results for the outer solution in the interior and exterior
domains (see Eq. (6.11) of Part I) as we approach the mouths in the form

$\displaystyle\lim_{\theta\rightarrow\theta_{0}}\lim_{r\rightarrow
b}\phi_{\mathrm{ext}}$
$\displaystyle\sim\dfrac{2\mathrm{i}A}{\pi}\left[\gamma_{\mathrm{e}}-\dfrac{\mathrm{i}\pi}{2}+\log\left(\dfrac{\widetilde{r}}{2}\right)\right]+\sum\limits_{n=-\infty}^{\infty}b_{n}\mathrm{Y}_{n}(b)\mathrm{e}^{\mathrm{i}n\theta_{0}}$
$\displaystyle\hskip
85.35826pt-\sum\limits_{n=-\infty}^{\infty}\left\\{\dfrac{AQ_{n}}{2}+b_{n}\mathrm{Y}_{n}^{\prime}(b)\mathrm{e}^{\mathrm{i}n\theta_{0}}\right\\}\dfrac{\mathrm{J}_{n}(b)}{\mathrm{J}_{n}^{\prime}(b)},$
(9a) $\displaystyle\lim_{\theta\rightarrow\theta_{0}}\lim_{r\rightarrow
a}\phi_{\mathrm{int}}$
$\displaystyle\sim\dfrac{2\mathrm{i}B}{\pi}\left[\gamma_{\mathrm{e}}-\dfrac{\mathrm{i}\pi}{2}+\log\left(\dfrac{\check{r}}{2}\right)\right]-\dfrac{B}{2}\sum\limits_{n=-\infty}^{\infty}\dfrac{\check{Q}_{n}}{\mathrm{J}_{n}^{\prime}(a)}\mathrm{J}_{n}(a),$
(9b)

where
$Q_{n}=\mathrm{J}_{n}(b)\mathrm{H}_{n}^{(1)\prime}(b)+\mathrm{J}_{n}^{\prime}(b)\mathrm{H}_{n}^{(1)}(b)$
and
$\check{Q}_{n}=\mathrm{J}_{n}(a)\mathrm{H}_{n}^{(1)\prime}(a)+\mathrm{J}_{n}^{\prime}(a)\mathrm{H}_{n}^{(1)}(a)$,
along with
$\widetilde{r}=\sqrt{(x-b\cos\theta_{0})^{2}+(y-b\sin\theta_{0})^{2}}$,
$\check{r}=\sqrt{(x-a\cos\theta_{0})^{2}+(y-a\sin\theta_{0})^{2}}$, and
$a=b-2\mathbcal{m}$. In the above, $\mathrm{J}_{n}(z)$ and $\mathrm{Y}_{n}(z)$
are Bessel functions of the first and second kind, respectively, and
$\mathrm{H}_{n}^{(1)}(z)$ are Hankel functions of the first kind. Having
obtained asymptotic forms for the outer solution in the interior, neck, and
exterior regions, we now consider the task of determining inner solutions at
the junctions to the resonator neck.

###  2.3 Inner solutions and asymptotic matching procedure

As before, we first rotate and translate the lattice so that the exterior
mouth of the resonator is located at the origin in
$(\widetilde{x},\widetilde{y})$ coordinates. We now introduce the inner
scaling $\widetilde{X}=\widetilde{x}/\varepsilon$ and
$\widetilde{Y}=\widetilde{y}/\varepsilon$ along with a regular expansion for
$\phi$ as in Part I; substituting these into the system we obtain the leading-
order inner boundary value problem

$\displaystyle(\partial_{\widetilde{X}}^{2}+\partial_{\widetilde{Y}}^{2})\Phi=0,\quad\mbox{
for }\widetilde{X}\in\widetilde{S}_{\mathrm{M}},$ (10a)
$\displaystyle\partial_{\widetilde{N}}\Phi=0,\quad\mbox{ for
}\widetilde{X}\in\partial\widetilde{S}_{\mathrm{M}},$ (10b)

where $\partial_{\widetilde{N}}$ denotes the normal derivative,
$\widetilde{S}_{\mathrm{M}}=\\{(\widetilde{X},\widetilde{Y}):-\infty<\widetilde{Y}<0,|\widetilde{X}|\leq
1\\}\cup\\{(\widetilde{X},\widetilde{Y}):\widetilde{Y}>0\\}$ is the exterior
mouth domain shown in Figure 2, and for clarity we omit the subscript for
$\Phi_{0}$, which denotes the leading term in the inner expansion as
$\varepsilon\rightarrow 0$. To obtain a solution we map the exterior mouth
region $\widetilde{S}_{\mathrm{M}}$ to the upper-half plane via the
Schwartz–Christoffel mapping

$\widetilde{Z}(W)=\frac{2}{\pi}\left[(W^{2}-1)^{1/2}-\mathrm{i}\log\left((W^{2}-1)^{1/2}+\mathrm{i}\right)+\mathrm{i}\log
W\right],$ (11)

where $\widetilde{Z}=\widetilde{X}+\mathrm{i}\widetilde{Y}$ and
$W=U+\mathrm{i}V$, which exhibits the asymptotic behaviours

$\displaystyle\lim_{W\rightarrow 0}\widetilde{Z}(W)\sim
1+\frac{2\mathrm{i}}{\pi}(1-\log 2)+\frac{2\mathrm{i}}{\pi}\log W,\quad\mbox{
and }\quad\lim_{W\rightarrow\infty}\widetilde{Z}(W)\sim\frac{2W}{\pi}.$ (12a)

The appropriate solution to Laplace’s equation in the upper-half plane
satisfying Neumann conditions along $V=0$ is given by

$\Phi=C_{1}\mathrm{Re}(\log W)+C_{2},$ (13)

where $C_{1}$ and $C_{2}$ are as yet unknown and from (12) it follows that

$\displaystyle\lim_{\begin{subarray}{c}\widetilde{Z}\rightarrow\infty\\\
(UHP)\end{subarray}}\Phi=C_{1}\log\left(\frac{\pi\widetilde{R}}{2}\right)+C_{2},\quad\mbox{
and }\quad\lim_{\begin{subarray}{c}\widetilde{Y}\rightarrow-\infty\\\
|\widetilde{X}|<1\end{subarray}}\Phi=C_{1}\left[\frac{\pi\widetilde{Y}}{2}-1+\log
2\right]+C_{2}.$ (14a) For the interior mouth region, the inner solution
$\Psi$ is obtained by a treatment analogous to that outlined above, but now
expressed in terms of the inner transformation
$\check{X}=\check{x}/\varepsilon$ and
$\check{Y}=\check{y}/\varepsilon=(\widetilde{y}+2\mathbcal{m})/\varepsilon$.
Accordingly, in this region, from (11), we expect the asymptotic behaviour
$\displaystyle\lim_{\begin{subarray}{c}\check{Z}\rightarrow\infty\\\
(LHP)\end{subarray}}\Psi=C_{3}\log\left(\frac{\pi\check{R}}{2}\right)+C_{4},\quad\mbox{
and }\quad\lim_{\begin{subarray}{c}\check{Y}\rightarrow\infty\\\
|\check{X}|<1\end{subarray}}\Psi=C_{3}\left[-\frac{\pi\check{Y}}{2}-1+\log
2\right]+C_{4},$ (14b)

where
$\check{Z}=\check{X}+\mathrm{i}\check{Y}=\check{R}\,\mathrm{exp}(\mathrm{i}\check{\Theta})$
and both $C_{3}$ and $C_{4}$ are as yet unknown. Subsequently, the matching
procedure [11], at leading order, gives rise to the relations

$\displaystyle\lim_{\begin{subarray}{c}\widetilde{R}\rightarrow\infty\\\
(UHP)\end{subarray}}\Phi\bigg{|}_{\widetilde{R}=\widetilde{r}/\varepsilon}=\lim_{\widetilde{r}\rightarrow
0}\lim_{\theta\rightarrow\theta_{0}}\phi_{\mathrm{ext}},\quad$
$\displaystyle\mbox{and}\quad\lim_{\widetilde{Y}\rightarrow-\infty}\Phi\bigg{|}_{\widetilde{X}=\widetilde{x}/\varepsilon,\widetilde{Y}=\widetilde{y}/\varepsilon}\qquad=\lim_{\widetilde{x}\rightarrow
0}\lim_{\widetilde{y}\rightarrow 0}\phi_{\mathrm{neck}},$ (15a)
$\displaystyle\lim_{\begin{subarray}{c}\check{R}\rightarrow\infty\\\
(LHP)\end{subarray}}\Psi\bigg{|}_{\check{R}=\check{r}/\varepsilon}=\lim_{\check{r}\rightarrow
0}\lim_{\theta\rightarrow\theta_{0}}\phi_{\mathrm{int}},\quad$
$\displaystyle\mbox{and}\quad\lim_{\check{Y}\rightarrow\infty}\Psi\bigg{|}_{\check{X}=\check{x}/\varepsilon,\check{Y}=\check{y}/\varepsilon}=\lim_{\check{x}\rightarrow
0}\lim_{\check{y}\rightarrow 0}\phi_{\mathrm{neck}}.$ (15b)

Thus, matching polynomial orders between (8) and (14), in addition to
logarithmic and non-logarithmic terms between (9) and (14), allows us to
determine, after significant algebra, all coefficients $C_{j}$, as well as the
monopole amplitudes
$B=2\mathrm{i}A\,(2\mathrm{i}\tau_{1}/\pi+\tau_{2}\tau_{5})^{-1}/\pi$ and

$A=\frac{2}{\pi
b\bar{\bar{h}}_{\varepsilon}}\sum\limits_{n=-\infty}^{\infty}g_{n},$ (16)

where
$g_{n}=b_{n}\mathrm{e}^{\mathrm{i}n\theta_{0}}/\mathrm{J}_{n}^{\prime}(b)$,
and

$\bar{\bar{h}}_{\varepsilon}=\dfrac{2\mathrm{i}}{\pi}\left[\gamma_{\mathrm{e}}-\dfrac{\mathrm{i}\pi}{2}-\log\left(\frac{\pi}{\varepsilon}\right)-\left(\frac{2\mathrm{i}}{\pi}\tau_{3}+\tau_{4}\tau_{5}\right)\left(\frac{2\mathrm{i}}{\pi}\tau_{1}+\tau_{2}\tau_{5}\right)^{-1}\right]-\frac{1}{2}\sum\limits_{n=-\infty}^{\infty}\dfrac{Q_{n}\mathrm{J}_{n}(b)}{\mathrm{J}_{n}^{\prime}(b)},$
(17)

along with

$\displaystyle\tau_{1}$ $\displaystyle=\frac{2\varepsilon}{\pi}(1-\log
2)\sin(2\mathbcal{m})-\cos(2\mathbcal{m}),\qquad\tau_{4}=-\frac{2\varepsilon}{\pi}(1-\log
2)\sin(2\mathbcal{m})+\cos(2\mathbcal{m}),$ (18a) $\displaystyle\tau_{2}$
$\displaystyle=-\frac{2\varepsilon}{\pi}\sin(2\mathbcal{m}),\hskip
56.9055pt\tau_{5}=\dfrac{2\mathrm{i}}{\pi}\left[\gamma_{\mathrm{e}}-\dfrac{\mathrm{i}\pi}{2}-\log\left(\frac{\pi}{\varepsilon}\right)\right]-\dfrac{1}{2}\sum\limits_{n=-\infty}^{\infty}\dfrac{\check{Q}_{n}}{\mathrm{J}_{n}^{\prime}(a)}\mathrm{J}_{n}(a),$
(18b) $\displaystyle\tau_{3}$
$\displaystyle=\left[\frac{2\varepsilon}{\pi}(1-\log
2)^{2}-\frac{\pi}{2\varepsilon}\right]\sin(2\mathbcal{m})-2(1-\log
2)\cos(2\mathbcal{m}).$ (18c)

Hence we obtain the same multipole eigensystem as in Eq. (4.16) of Part I,
which we express as

$\dfrac{\mathrm{i}A}{2}\left(\frac{\mathrm{J}_{n}^{\prime}(b)\mathrm{Y}_{n}(b)+\mathrm{Y}_{n}^{\prime}(b)\mathrm{J}_{n}(b)}{\mathrm{J}_{n}^{\prime}(b)\mathrm{Y}_{n}^{\prime}(b)}+2\sum_{m=-\infty}^{\infty}(-1)^{n+m}S_{m-n}^{\mathrm{Y}}(\mathbf{k}_{\mathrm{B}})\frac{\mathrm{J}_{m}(b)}{\mathrm{Y}_{n}^{\prime}(b)}\mathrm{e}^{-\mathrm{i}(m-n)\theta_{0}}\right)\\\
+g_{n}+\sum_{m=-\infty}^{\infty}(-1)^{m+n}S_{m-n}^{\mathrm{Y}}(\mathbf{k}_{\mathrm{B}})\frac{\mathrm{J}_{m}^{\prime}(b)}{\mathrm{Y}_{n}^{\prime}(b)}\mathrm{e}^{-\mathrm{i}(m-n)\theta_{0}}g_{m}=0,$
(19)

but with the replacement $h_{\varepsilon}\mapsto\bar{\bar{h}}_{\varepsilon}$
as written in (16) above. Thus, it would appear that we have an eigensystem
for examining specially-scaled resonator arrays, however, it turns out that
the above formulation exhibits pathological behaviour at frequencies near
$\bar{\bar{h}}_{\varepsilon}=0$ which occurs at much lower frequencies than
for the wall thickness values discussed in Part I. We now present a
regularisation procedure for resolving this issue.

##  3 Regularised multipole system formulation

As indicated in the previous section, although the system in Part I with the
replacement $h_{\varepsilon}\mapsto\bar{\bar{h}}_{\varepsilon}$ is indeed
formally correct, it can return spurious spectral behaviours upon truncating
for numerical evaluation (i.e., we may observe incorrect folded bands, and
flat band surfaces, that are not part of the genuine spectrum in a
neighbourhood around frequencies corresponding to
$\bar{\bar{h}}_{\varepsilon}=0$). We easily resolve this issue by summing (19)
over all $n$ to obtain a more numerically stable yet formally equivalent
representation for $A$ (c.f., to the form (16) above) as

$A=2\mathrm{i}\left(E+2\left[\sum_{v=-\infty}^{\infty}J_{v}(b)F_{v}\right]-\mathrm{i}\pi
b\bar{\bar{h}}_{\varepsilon}\right)^{-1}\sum_{m=-\infty}^{\infty}\mathrm{J}_{m}^{\prime}(b)F_{m}g_{m},$
(20)

where

$E=\sum_{u=-\infty}^{\infty}E_{u},\quad\mbox{and}\quad
F_{m}=\sum_{p=-\infty}^{\infty}\dfrac{(-1)^{p-m}S_{m-p}^{\mathrm{Y}}(\mathbf{k}_{\mathrm{B}})\mathrm{e}^{\mathrm{i}(p-m)\theta_{0}}}{\mathrm{Y}_{p}^{\prime}(b)},$
(21)

with
$E_{u}=\left[\mathrm{J}_{u}^{\prime}(b)\mathrm{Y}_{u}(b)+\mathrm{Y}_{u}^{\prime}(b)\mathrm{J}_{u}(b)\right]/\left[\mathrm{J}_{u}^{\prime}(b)\mathrm{Y}_{u}^{\prime}(b)\right]$.
Substituting this new representation for $A$ (20) into (19) admits the
regularised system

$g_{n}+\sum_{m=-\infty}^{\infty}(-1)^{m+n}S_{m-n}^{\mathrm{Y}}(\mathbf{k}_{\mathrm{B}})\frac{\mathrm{J}_{m}^{\prime}(b)}{\mathrm{Y}_{n}^{\prime}(b)}\mathrm{e}^{-\mathrm{i}(m-n)\theta_{0}}g_{m}-\dfrac{\chi_{n}}{H_{\varepsilon}}\sum_{m=-\infty}^{\infty}\mathrm{J}_{m}^{\prime}(b)F_{m}g_{m}=0,$
(22)

for all $n$, where

$\chi_{n}=E_{n}+2\sum_{p=-\infty}^{\infty}(-1)^{n+p}S_{p-n}^{\mathrm{Y}}(\mathbf{k}_{\mathrm{B}})\frac{\mathrm{J}_{p}(b)}{\mathrm{Y}_{n}^{\prime}(b)}\mathrm{e}^{-\mathrm{i}(p-n)\theta_{0}},$
(23a) and
$H_{\varepsilon}=E+2\left[\sum_{v=-\infty}^{\infty}J_{v}(b)F_{v}\right]-\mathrm{i}\pi
b\bar{\bar{h}}_{\varepsilon}.$ (23b)

On examining the system (22) numerically, we find that spurious effects are
removed and the genuine spectrum is observed, see for example Figure 3 which
is discussed in further detail below. In order for the regularisation to be
effective, we remark that the known expressions $E$, $\chi_{n}$, $F_{m}$, and
$H_{\varepsilon}$ must be suitably converged with appropriately chosen
truncation numbers for the sums, which we represent by $L_{E}$,
$L_{\chi_{n}}$, $L_{F_{m}}$, and $L_{H_{\varepsilon}}$.

##  4 Asymptotic representations for the dispersion equation

In this section, we construct an asymptotic representation of the dispersion
equation for the regularised system (22) at low frequencies. We begin by
observing that the lattice sums $S_{m}^{\mathrm{Y}}$, which feature in (22) as
well as in the functions $F_{m}$ (21), $\chi_{n}$ (23a), and $H_{\varepsilon}$
(23b), can be expressed as

$S_{m}^{\mathrm{Y}}\approx\sum_{r=-\Omega_{m}}^{\infty}\beta_{r}^{(m)}b^{r},$
(24)

where
$\Omega_{m}=2\left[\mathrm{Floor}\left\\{(m-3)/4\right\\}+\mathrm{Floor}\left\\{(m-4)/4\right\\}+2\right]$
for $m\geq 4$, denotes the first non-zero order of the lattice sum (i.e.,
$\Omega_{4}=4$ and $\Omega_{7}=6$), along with
$\Omega_{0}=\Omega_{1}=\Omega_{2}=\Omega_{3}=2$ [13, A004524]. Closed-form
expressions for the first few lattice sums, i.e., the values $\beta_{r}^{(m)}$
for small $m$, are derived and presented in Appendix A for reference. Note
that these terms may contain logarithmic behaviour in the lattice spacing $d$
and therefore in $b$ if the area fraction $f=\pi b^{2}/d^{2}$ is fixed. With
the lattice sum expansions (24), we examine the system (22) in the limit of
vanishing $b$ and to within a dipole truncation in $n$, which admits the
matrix form

$(H_{\varepsilon}\mathbf{C}-\mathbf{P}-\mathbf{Q})\mathbf{g}\approx\boldsymbol{0},$
(25)

where $\mathbf{C}$ is the matrix representation of the system for an array of
Neumann cylinders (obtained by setting $\chi_{n}=0$ in (22), see [14, Eq.
(3.120)]), $\mathbf{P}$ and $\mathbf{Q}$ are perturbation matrices due to the
(specially-scaled) aperture, and $\mathbf{g}$ is the vector of $g_{m}$
coefficients. These matrices take the forms

$\mathbf{C}\approx\left[\begin{array}[]{c|c|c}\phantom{\bigg{|}^{1}}1+\dfrac{\pi}{4}\beta_{-2}^{(0)}&\dfrac{\pi
b}{4}\mathrm{e}^{\mathrm{i}\theta_{0}}\left(\beta_{-2}^{(1)}\right)^{\ast}&-\dfrac{\pi}{4}\mathrm{e}^{2\mathrm{i}\theta_{0}}\left(\beta_{-2}^{(2)}\right)^{*}\\\
\hline\cr\phantom{\bigg{|}^{1}}-\dfrac{\pi}{4b}\mathrm{e}^{-\mathrm{i}\theta_{0}}\beta_{-2}^{(1)}&1-\dfrac{\pi}{4}\beta_{-2}^{(0)}&\dfrac{\pi}{4b}\mathrm{e}^{\mathrm{i}\theta_{0}}\left(\beta_{-2}^{(1)}\right)^{*}\\\
\hline\cr-\dfrac{\pi}{4}\mathrm{e}^{-2\mathrm{i}\theta_{0}}\beta_{-2}^{(2)}&-\dfrac{\pi
b}{4}\mathrm{e}^{-\mathrm{i}\theta_{0}}\beta_{-2}^{(1)}&1+\dfrac{\pi}{4}\beta_{-2}^{(0)}\end{array}\right],$
(26a)
$\mathbf{P}\approx\left[\begin{array}[]{c|c|c}\phantom{\bigg{|}^{1}}\dfrac{\pi^{2}\left|\beta_{-2}^{(1)}\right|^{2}}{4b}&\dfrac{\pi^{2}\mathrm{e}^{\mathrm{i}\theta_{0}}(\beta_{-2}^{(1)})^{\ast}\beta_{-2}^{(0)}}{4}&-\dfrac{\pi^{2}\mathrm{e}^{2\mathrm{i}\theta_{0}}((\beta_{-2}^{(1)})^{\ast})^{2}}{4b}\\\
\hline\cr\phantom{\bigg{|}^{1}}-\dfrac{\pi\left(\pi\beta_{-2}^{(0)}-2\right)\mathrm{e}^{-\mathrm{i}\theta_{0}}\beta_{-2}^{(1)}}{4b^{2}}&-\dfrac{\pi\beta_{-2}^{(0)}\left(\pi\beta_{-2}^{(0)}-2\right)}{4b}&\dfrac{\pi\left(\pi\beta_{-2}^{(0)}-2\right)\mathrm{e}^{\mathrm{i}\theta_{0}}(\beta_{-2}^{(1)})^{\ast}}{4b^{2}}\\\
\hline\cr\phantom{\bigg{|}^{1}}-\dfrac{\pi^{2}\mathrm{e}^{-2\mathrm{i}\theta_{0}}(\beta_{-2}^{(1)})^{2}}{4b}&-\dfrac{\pi^{2}\mathrm{e}^{-\mathrm{i}\theta_{0}}\beta_{-2}^{(0)}\beta_{-2}^{(1)}}{4}&\dfrac{\pi^{2}\left|\beta_{-2}^{(1)}\right|^{2}}{4b}\end{array}\right],$
(26b) and
$\mathbf{Q}\approx\left[\begin{array}[]{c|c|c}\phantom{\bigg{|}^{1}}\dfrac{\pi^{2}\mathrm{e}^{-\mathrm{i}\theta_{0}}\beta_{-2}^{(1)}}{2}\eta^{\ast}+\mathrm{e}^{\mathrm{i}\theta_{0}}\pi^{2}(\beta_{-2}^{(1)})^{\ast}\eta&\dfrac{\pi^{2}b}{2}\beta_{-2}^{(0)}\eta^{\ast}&\dfrac{\mathrm{e}^{\mathrm{i}\theta_{0}}\pi^{2}(\beta_{-2}^{(1)})^{\ast}}{2}\eta^{\ast}\\\
\hline\cr\phantom{\bigg{|}^{1}}\dfrac{\pi\eta}{b}\left(2-\pi\beta_{-2}^{(0)}\right)&0&\dfrac{\pi\eta^{\ast}}{b}\left(2-\pi\beta_{-2}^{(0)}\right)\\\
\hline\cr-\dfrac{\pi^{2}\mathrm{e}^{-\mathrm{i}\theta_{0}}\beta_{-2}^{(1)}}{2}\eta&\dfrac{\pi^{2}b}{2}\beta_{-2}^{(0)}\eta&-\dfrac{\pi^{2}\mathrm{e}^{\mathrm{i}\theta_{0}}(\beta_{-2}^{(1)})^{\ast}}{2}\eta-\mathrm{e}^{-\mathrm{i}\theta_{0}}\pi^{2}\beta_{-2}^{(1)}\eta^{\ast}\end{array}\right],$
(26c)

where we remark that next order terms are taken in most matrix entries above
for our analysis in the following section, but are omitted here for
compactness, and

$\eta=-\frac{\beta_{-2}^{(0)}}{4}+\frac{1}{4}\mathrm{e}^{-2\mathrm{i}\theta_{0}}\beta_{-2}^{(2)}+\sum_{q=1}^{\infty}\left\\{\frac{\beta_{-2q-2}^{(2q+2)}\mathrm{e}^{-(2q+2)\mathrm{i}\theta_{0}}}{4^{q+1}(2q+1)!}\right\\},$
(27)

with $\ast$ representing the complex conjugate operation. In addition, we have

$H_{\varepsilon}\approx\frac{-2+4\bar{\bar{f}}_{\varepsilon}+\pi\beta_{-2}^{(0)}}{b}-\frac{\pi}{2}\left[\mathrm{e}^{\mathrm{i}\theta_{0}}(\beta_{-2}^{(1)})^{\ast}-\beta_{-2}^{(1)}\mathrm{e}^{-\mathrm{i}\theta_{0}}\right]+O(b),$
(28)

where we take an analogous scaling to that found in Part I

$\bar{\bar{f}}_{\varepsilon}=\frac{\pi
b^{2}}{4\mathrm{i}}\bar{\bar{h}}_{\varepsilon},$ (29a) as well as a dipole
truncation in $\bar{\bar{h}}_{\varepsilon}$ (17) to obtain the asymptotic form
$\lim_{b\rightarrow 0}\lim_{a\rightarrow 0}\lim_{\mathbcal{m}\rightarrow
0}\bar{\bar{h}}_{\varepsilon}\sim\dfrac{2\mathrm{i}}{\pi}\left\\{\frac{1}{b^{2}}-\frac{1}{8}-\log\left(\frac{\pi
b}{2\varepsilon}\right)\right.\\\
\left.+\left[\frac{1}{a^{2}}-\frac{\mathbcal{m}\pi}{\varepsilon}-\frac{17}{8}-\log\left(\frac{\pi
a}{8\varepsilon}\right)\right]\left[1+\frac{4\varepsilon\mathbcal{m}}{\pi}\left(\frac{1}{a^{2}}-\frac{1}{8}-\log\left(\frac{\pi
a}{2\varepsilon}\right)\right)\right]^{-1}\right\\}.$ (29b)

With the above expansions, we are now able to construct closed-form
representations for the dispersion equation of specially-scaled resonators,
however before proceeding to this task, we comment that certain
$\beta_{-2n}^{(2n)}$ terms in (27) are vanishing for square lattices, such as
$\beta_{-6}^{(6)}$, $\beta_{-10}^{(10)}$, and $\beta_{-14}^{(14)}$ and that
many terms are readily extracted from the explicit forms in Appendix A.

###  4.1 Dispersion equation forms along selected symmetry directions

It is instructive to consider cases where the resonator geometry has natural
symmetry. Thus, with the asymptotic system (25) in mind, we now take
$\theta_{0}=0$ and consider Bloch coordinates located on the high symmetry
planes of the Brillouin zone, where $(k_{\mathrm{B}},\theta_{\mathrm{B}})$
denotes the polar form of the dimensionless Bloch vector
$\mathbf{k}_{\mathrm{B}}$.

For reference, an outline of the symmetry planes for the fundamental cell in
reciprocal space (Brillouin zone) is presented in Figure 1 of Part I, showing
the paths between the high symmetry Bloch vector coordinates $\Gamma=(0,0)$,
$X=(\pi/d,0)$, $M=(\pi/d,\pi/d)$, and $Y=(0,\pi/d)$. In the first instance, we
evaluate the determinant of the system (25) along the $\Gamma X$ (i.e.,
$\theta_{\mathrm{B}}=0$) direction, which to leading order, returns the
isotropic result from Part I (Eq. (4.22)) with the replacement
$f_{\varepsilon}\mapsto\bar{\bar{f}}_{\varepsilon}$ in the form

$k_{\mathrm{B}}^{2}=\frac{1+f}{1-f}\left(1-\frac{2f(1-\bar{\bar{f}}_{\varepsilon})}{1-2\bar{\bar{f}}_{\varepsilon}}\right),$
(30)

where $f=\pi b^{2}/d^{2}$ denotes the area ratio and
$\bar{\bar{f}}_{\varepsilon}$ is given in (29a). In contrast, along the
$\Gamma Y$ direction ($\theta_{\mathrm{B}}=\pi/2$) we observe a much more
complicated form, which emerges from the symmetries of the system matrices

$\mathbf{C}=\left[\begin{array}[]{rrr}\bullet&\bowtie&\boxtimes\\\
\diamond&\circledcirc&-\diamond\\\
\boxtimes&-\bowtie&\bullet\end{array}\right],\quad\mathbf{P}=\left[\begin{array}[]{rrr}\bullet&\bowtie&-\bullet\\\
\diamond&\circledcirc&-\diamond\\\
-\bullet&-\bowtie&\bullet\end{array}\right],\quad\mathbf{Q}=\left[\begin{array}[]{ccc}\bullet&\bowtie&-\boxtimes\\\
\diamond&0&\diamond\\\ \boxtimes&\bowtie&-\bullet\end{array}\right],$

where symbols are used to represent the symmetry of each matrix independently
(and does not imply equivalence of values between matrices). Subsequently,
after a significant amount of algebra, we obtain the dispersion equation

$\left(\sqrt{\frac{a^{2}\left(2f+(1-2\bar{\bar{f}}_{\varepsilon})k_{\mathrm{B}}^{2}\right)+b^{2}\left(k_{\mathrm{B}}^{2}+1\right)-\Lambda_{\Gamma\mathrm{Y}}}{2a^{2}(1-2\bar{\bar{f}}_{\varepsilon})+b^{2}}}-1\right)\cdot\\\
\left(\sqrt{\frac{a^{2}\left(f^{2}+f(2\bar{\bar{f}}_{\varepsilon}-1)k_{\mathrm{B}}^{2}+f+(1-2\bar{\bar{f}}_{\varepsilon})k_{\mathrm{B}}^{2}\right)+b^{2}\left(-f^{2}-fk_{\mathrm{B}}^{2}+k_{\mathrm{B}}^{2}+1\right)+\xi_{\Gamma\mathrm{Y}}}{2\left(1-f^{2}\right)\left(a^{2}(1-2\bar{\bar{f}}_{\varepsilon})+b^{2}\right)}}-1\right)=0,$
(31) where
$\Lambda_{\Gamma\mathrm{Y}}^{2}=a^{4}(2f+(1-2\bar{\bar{f}}_{\varepsilon})k_{\mathrm{B}}^{2})^{2}+2a^{2}b^{2}[2f(k_{\mathrm{B}}^{2}+1)+(1-2\bar{\bar{f}}_{\varepsilon})\left(k_{\mathrm{B}}^{2}-1\right)k_{\mathrm{B}}^{2}]+b^{4}(k_{\mathrm{B}}^{2}-1)^{2},$
(32a) and
$\xi_{\Gamma\mathrm{Y}}^{2}=\left[a^{2}\left(f^{2}+f(2\bar{\bar{f}}_{\varepsilon}-1)k_{\mathrm{B}}^{2}+f+(1-2\bar{\bar{f}}_{\varepsilon})k_{\mathrm{B}}^{2}\right)+b^{2}\left(-f^{2}-fk_{\mathrm{B}}^{2}+k_{\mathrm{B}}^{2}+1\right)\right]^{2}\\\
-4b^{2}(f-1)\left(f^{2}-1\right)k_{\mathrm{B}}^{2}\left(a^{2}(1-2\bar{\bar{f}}_{\varepsilon})+b^{2}\right).$
(32b)

The key point here is that although the dispersion equation along
$\Gamma\mathrm{Y}$ in (31) above is indeed accurate, it is however intractable
and does not provide any useful insight, despite its correspondence to the
high symmetry settings of $\theta_{0}=0$ and $\theta_{\mathrm{B}}=\pi/2$.
Attempts in reducing (31), for example, by taking a Taylor series expansion in
small $k_{\mathrm{B}}$, provide equally unhelpful forms. Despite its length,
the expression (31) does however demonstrate that incorporating anisotropy to
the dispersion equation for Helmholtz resonator arrays in the specially-scaled
limit is a much more formidable task than for the (moderately) thin-walled
case, as shown in Part I.

The unexpected complexity of the dispersion equation (31) can be understood by
considering the geometry of the resonator: in the specially-scaled limit, the
neck of each resonator is so thin and long that the resonator is almost
invariant under all rotation and reflection operations for the square lattice.
Subsequently the metamaterial may be considered almost-isotropic, even as we
approach the resonance frequency, and so next-order asymptotic corrections
take on much more complicated forms. Numerical investigations confirm much
smaller anisotropy in this scaled setting, to the point where the medium may
be essentially regarded as isotropic for practical applications. Accordingly,
the low-frequency dispersion equation for these specially-scaled resonators
may be taken to be the isotropic form (30) (which in fact holds for all values
of $\theta_{0}$) over the entirety of the Brillouin zone. The effectiveness of
this approximation is examined numerically in Section 7 below.

##  5 Helmholtz resonance (cut-off) condition

As determined in Part I and from the denominator of (30) above, the Helmholtz
resonance condition in the specially-scaled limit is given by
$1-2\bar{\bar{f}}_{\varepsilon}\approx 0$ and requires careful examination, as
such resonators are capable of achieving very low-frequency resonances.
Accordingly, we return to the asymptotic form of $\bar{\bar{h}}_{\varepsilon}$
in (29b) and subsequently write $1-2\bar{\bar{f}}_{\varepsilon}\approx 0$ in
the form

$\frac{1}{8}+\log\left(\frac{\pi
b}{2\varepsilon}\right)-\left[\frac{1}{a^{2}}-\frac{\mathbcal{m}\pi}{\varepsilon}-\frac{17}{8}-\log\left(\frac{\pi
a}{8\varepsilon}\right)\right]\left[1+\frac{4\varepsilon\mathbcal{m}}{\pi}\left(\frac{1}{a^{2}}\right)\right]^{-1}\approx
0,$ (33)

where we remark that $\tau_{5}$ in (18) is taken to $O(1)$ in the numerator
and to leading order in the denominator. For the purposes of analysis, we do
not advise solving (33) above to determine the conditions for resonance, as it
is unclear which terms play a leading role in the small $\varepsilon$ limit;
we seek the preferred scalings between $\varepsilon$, $a$ (or $b$), and
$\mathbcal{m}$ that give the lowest frequency resonance.

Therefore, it is helpful to introduce the scalings
$\mathbcal{m}=\kappa_{\mathbcal}{m}\varepsilon^{\mu}$ (i.e., where $0<\mu<1$,
$\varepsilon\rightarrow 0$, and $\kappa_{\mathbcal}{m}=O(1)$), and
$a=\kappa_{a}\varepsilon^{\gamma}$ (i.e., where $0<\gamma<1$,
$\varepsilon\rightarrow 0$, and $\kappa_{a}=O(1)$) which admits

$\frac{1}{8}+\log\left(\frac{\pi\left[\kappa_{a}\varepsilon^{\gamma-1}+2\kappa_{\mathbcal}{m}\varepsilon^{\mu-1}\right]}{2}\right)-\dfrac{\left[\dfrac{1}{\kappa_{a}^{2}}\varepsilon^{-2\gamma}-\kappa_{\mathbcal}{m}\varepsilon^{\mu-1}\pi-\dfrac{17}{8}-\log\left(\dfrac{\pi\kappa_{a}\varepsilon^{\gamma-1}}{8}\right)\right]}{\left[1+\dfrac{4\kappa_{\mathbcal}{m}}{\pi\kappa_{a}^{2}}\varepsilon^{1+\mu-2\gamma}\right]}\approx
0.$ (34)

The radii $a=\kappa_{a}\varepsilon^{\gamma}$ and
$b=\kappa_{a}\varepsilon^{\gamma}+2\kappa_{\mathbcal}{m}\varepsilon^{\mu}$,
must be of the same order, requiring $\mu\geq\gamma$. We now examine, in some
detail, the different dominant balance scalings that are possible for the
representation (34).

####  5.0.1 Dominant balance in all three numerator terms

From the form (34), we take the scaling $1+\mu-2\gamma>0$, where to ensure
that the $O(\varepsilon^{-2\gamma})$, $O(\varepsilon^{\mu-1})$, and $O(1)$
terms balance in the numerator, we require $\gamma=0$ and $\mu=1$.
Accordingly, the Helmholtz resonance condition then takes the form

$\frac{9}{4}+2\log\left(\frac{\pi\kappa_{a}}{4\varepsilon}\right)\approx\left[\frac{1}{\kappa_{a}^{2}}-\pi\kappa_{\mathbcal}{m}\right],$
(35a) as $\varepsilon\to 0$. Expressing the above in dimensional terms we find
$k_{\mathrm{max}}\approx\frac{1}{\bar{a}}\sqrt{\dfrac{1}{9/4+2\log\left(\pi\bar{a}/(4\bar{\ell})\right)+\pi\bar{\mathbcal{m}}/\bar{\ell}}},$
(35b)

as the Helmholtz resonance condition under this dominant balance scaling.

####  5.0.2 Dominant balance in numerator pairs

Another relationship emerges by considering (34) in the $1+\mu-2\gamma>0$
regime again, which gives

$\frac{9}{4}+\log\left(\frac{\pi^{2}\kappa_{a}\varepsilon^{\gamma-1}\left[\kappa_{a}\varepsilon^{\gamma-1}+2\kappa_{\mathbcal}{m}\varepsilon^{\mu-1}\right]}{16}\right)\approx\left[\frac{1}{\kappa_{a}^{2}}\varepsilon^{-2\gamma}-\pi\kappa_{\mathbcal}{m}\varepsilon^{\mu-1}\right],$
(36)

where balance on the right-hand side is achieved with the scaling
$\gamma=\frac{1}{2}(1-\mu)$, which, as $\mu\geq\gamma$, means that $\mu\geq
1/3$. Turning to the logarithmic argument we see that for $\mu>1/3$, the
resonance condition (36) takes the form

$\frac{9}{4}+\log\left(\frac{\pi^{2}\kappa_{a}^{2}\varepsilon^{-1-\mu}}{16}\right)\approx\left[\frac{1}{\kappa_{a}^{2}}-\pi\kappa_{\mathbcal}{m}\right]\varepsilon^{\mu-1},$
(37a) whereas for $\mu=1/3$ we observe
$\frac{9}{4}+\log\left(\frac{\pi^{2}\kappa_{a}\varepsilon^{-4/3}\left[\kappa_{a}+2\kappa_{\mathbcal}{m}\right]}{16}\right)\approx\left[\frac{1}{\kappa_{a}^{2}}-\pi\kappa_{\mathbcal}{m}\right]\varepsilon^{-2/3}.$
(37b)

Expressing these resonance conditions in terms of dimensional parameters we
again find for $\mu>1/3$ the form in (35b), and for $\mu=1/3$ the slightly
different form

$k_{\mathrm{max}}\approx\frac{1}{\bar{a}}\sqrt{\dfrac{1}{9/4+\log\left(\pi^{2}\bar{a}\bar{b}/(16\bar{\ell}^{2})\right)+\pi\bar{\mathbcal{m}}/\bar{\ell}}},$
(38)

with numerical investigations suggesting that only minor differences are found
between (35b) and (38), since $a$ and $b$ must be the same order. For
reference, we remark that under the dominant balance scaling
$\mu+1-2\gamma>0$, we may write

$\bar{\bar{f}}_{\varepsilon}\approx\frac{b^{2}}{2}\left[\frac{1}{a^{2}}+\frac{1}{b^{2}}-\frac{9}{4}-\frac{\pi\mathbcal{m}}{\varepsilon}-\log\left(\frac{\pi^{2}ab}{16\varepsilon^{2}}\right)\right].$
(39)

##  6 Homogenisation of Helmholtz resonator arrays

###  6.1 Classical homogenisation results for arrays of ideal cylinders

For an isotropic fluid medium of density $\rho$ and bulk modulus $B$,
structured with a two-dimensional array of isotropic fluid cylinders, of
density $\rho_{\mathrm{c}}$ and bulk modulus $B_{c}$, the effective density
and bulk modulus are given explicitly in the quasistatic limit by [15, 16]

$\rho_{\mathrm{eff}}^{-1}=\frac{1}{\rho}\frac{\rho_{c}(1-f)+\rho(1+f)}{\rho_{c}(1+f)+\rho(1-f)},\quad\mbox{and}\quad
B_{\mathrm{eff}}^{-1}=\frac{1-f}{B}+\frac{f}{B_{c}},$ (40)

where $f=\pi b^{2}/d^{2}$, as defined previously, denotes the filling
fraction. Subsequently, results corresponding to a two-dimensional array of
cylinders with Neumann boundary conditions on the walls take the form

$\displaystyle\rho_{\mathrm{eff}}^{-1}=\frac{1-f}{\rho(1+f)},\quad\mbox{and}\quad
B_{\mathrm{eff}}^{-1}=\frac{1-f}{B},$ (41a) which are obtained via the limits
$B_{c}\rightarrow\infty$ and $\rho_{c}\to\infty$ in (40) above. Substituting
(41a) into the (dimensional) dispersion equation for plane waves in an
unbounded isotropic medium
$\bar{k}_{\mathrm{B}}\,\rho_{\mathrm{eff}}^{-1}\,\bar{k}_{\mathrm{B}}-\omega^{2}B_{\mathrm{eff}}^{-1}=0,$
(41b) we obtain the (dimensionless) dispersion equation for an array of
Neumann cylinders $k_{\mathrm{B}}^{2}=1+f,$ (41c)

as seen from our analysis in Part I (Eq. (4.23)) and from Movchan et. al [14,
Eq. (3.158)], and where
$\bar{\mathbf{k}}_{\mathrm{B}}=k\mathbf{k}_{\mathrm{B}}$ is the nondimensional
Bloch vector. Similarly, if we return to the two-phase fluid array results
(40) and take the limits $B_{c}\rightarrow B$ and $\rho_{c}\to\infty$ we
obtain

$\displaystyle\rho_{\mathrm{eff}}^{-1}=\frac{1-f}{\rho(1+f)},\quad\mbox{and}\quad
B_{\mathrm{eff}}^{-1}=\frac{1}{B},$ (42a) which on substitution in (41b)
returns the dispersion relation $k_{\mathrm{B}}^{2}=\frac{1+f}{1-f},$ (42b)

and is identical to the lowest-order (isotropic) approximation for an array of
thin-walled Helmholtz resonators at quasi-static frequencies (i.e.,
$\omega\rightarrow 0$), as given in Part I and in Llewellyn–Smith [17]. That
is, the lowest-order (isotropic) approximation for the dispersion equation of
a Helmholtz resonator array (42b) at very low frequencies is indistinguishable
from an array of fluid cylinders that do not possess a contrast in bulk
modulus but are much denser than the background fluid. On comparing (41a) and
(42a) we see that the impact of introducing a small gap into the wall of a
perfect Neumann cylinder has the leading-order effect of modifying the
effective bulk modulus but not the effective density. As such, we may expect
to recover the same expressions for the bulk modulus and density in (42a) at
low frequencies from any isotropic descriptions for resonator arrays.

###  6.2 Isotropic descriptions for thin-walled Helmholtz resonator arrays

Having discussed the behaviour of the isotropic description at zero frequency,
we now restate the result obtained in Part I (Eq. (4.22))

$k_{\mathrm{B}}^{2}=\frac{1+f}{1-f}\left(1-\frac{2f(1-f_{\varepsilon})}{1-2f_{\varepsilon}}\right),$
(43)

which is valid for frequencies across the range of the first band surface,
where

$f_{\varepsilon}\sim 1-b^{2}/8+b^{2}\log\left(\theta_{\mathrm{ap}}/2\right),$
(44)

and $\theta_{\mathrm{ap}}$ denotes the aperture half-width angle. Here, we
emphasise that there are infinitely many ways in which this dispersion
equation may be decomposed into the form (41b). For example, we may extract
from (43) the response functions

$\rho_{\mathrm{eff}}^{-1}=\frac{(1-f)}{\rho(1+f)},\quad\mbox{and}\quad
B_{\mathrm{eff}}^{-1}=\frac{1}{B}\frac{(1-2f_{\varepsilon}-2f(1-f_{\varepsilon}))}{(1-2f_{\varepsilon})},$
(45)

as candidates for the behaviour of the resonator array, within an isotropic
approximation. Such expressions are consistent with the two-phase fluid
results (42a) as $\omega\rightarrow 0$, vanishing filling fraction results
$f\rightarrow 0$, and the expressions for closed cylinders (41a) as the
aperture is closed (i.e., as $f_{\varepsilon}\rightarrow\infty$). The absence
of dispersion in the density response is also consistent with existing
literature on the topic [1, 2]. Thus, we consider (45) as the effective
(homogenised) quantities of the medium for frequencies spanning the first band
surface, within an isotropic approximation.

###  6.3 Anisotropic descriptions for thin-walled Helmholtz resonator arrays

As established in Part I, thin-walled resonators generally exhibit strong
anisotropy at low frequencies, and so isotropic descriptions are insufficient
to accurately describe the first band surface. Accordingly, we restate the
anisotropic result for the first spectral band from Part I (Eq. (4.29a)):

$k_{\mathrm{B}}^{2}=\frac{(f+1)\left[b^{2}f(2f-1)+f_{\varepsilon}(f+1)(2f_{\varepsilon}f-2f_{\varepsilon}-2f+1)\right]}{b^{2}f\cos(2\left[\theta_{0}-\theta_{\mathrm{B}}\right])+b^{2}f^{2}+f_{\varepsilon}(2f_{\varepsilon}-1)\left(f^{2}-1\right)},$
(46)

and emphasise once more that the assignment of the effective density and bulk
modulus is non-unique in the dispersion equation for plane waves in an
anisotropic medium [18]

$\bar{k}_{\mathrm{B}i}\,(\rho_{\mathrm{eff}}^{-1})_{ij}\,\bar{k}_{\mathrm{B}j}-\omega^{2}B_{\mathrm{eff}}^{-1}=0.$
(47)

For example, we may propose the candidate forms

$\displaystyle\boldsymbol{\rho}^{-1}_{\mathrm{eff}}$
$\displaystyle=\frac{(1-f)}{\rho(1+f)}\left[\begin{array}[]{c|c}1-b^{2}f\left[1-\cos(2\theta_{0})\right]/\mathcal{H}&b^{2}f\sin(2\theta_{0})/\mathcal{H}\\\
\hline\cr
b^{2}f\sin(2\theta_{0})/\mathcal{H}&1-b^{2}f\left[1+\cos(2\theta_{0})\right]/\mathcal{H}\end{array}\right]$
(48c) $\displaystyle B_{\mathrm{eff}}^{-1}$
$\displaystyle=\frac{(1-f)}{B\mathcal{H}}\left[b^{2}f(2f-1)+(f+1)f_{\varepsilon}(2(f-1)f_{\varepsilon}-2f+1)\right],$
(48d)

where

$\mathcal{H}=b^{2}f^{2}+b^{2}f+(f^{2}-1)f_{\varepsilon}(2f_{\varepsilon}-1)=0.$
(49)

These expressions are consistent with both the isotropic results (45) and the
two-phase fluid results (42a) as $\omega\rightarrow 0$, as well as with the
closed aperture (Neumann cylinder) results (41a) as
$f_{\varepsilon}\rightarrow\infty$ and the vanishing fill fraction limit
$f\rightarrow 0$. Note that the expression $\mathcal{H}$ is the anisotropic
approximation to the Helmholtz resonance condition
$\det\\{(\rho_{\mathrm{eff}}^{-1})_{ij}\\}=0$.

###  6.4 Descriptions for moderately thick-walled Helmholtz resonator arrays

For moderately thick-walled resonators, the analysis proceeds as in Section
66.2 and 66.3 above but with the replacement
$f_{\varepsilon}\mapsto\check{f}_{\varepsilon}$, where the definition for
$\check{f}_{\varepsilon}$ is given in Part I (Eq. (6.14)).

###  6.5 Descriptions for specially-scaled Helmholtz resonator arrays

To obtain an isotropic description for the specially-scaled thick-walled
resonator configuration discussed in this Part II, the analysis proceeds as in
Section 66.2 with the replacement
$f_{\varepsilon}\mapsto\bar{\bar{f}}_{\varepsilon}$ where
$\bar{\bar{f}}_{\varepsilon}$ is defined in (39). As shown in Section 44.1
earlier, anisotropic descriptions derived from the regularised system for
specially-scaled resonators are generally intractable, i.e., see (31). That
said, it should be possible to obtain an anisotropic description of the first
spectral band alone using (46) with the replacement
$f_{\varepsilon}\mapsto\bar{\bar{f}}_{\varepsilon}$.

##  7 Numerical Results

In this section, we compute a broad selection of band diagrams, comparing
results from our regularised system (22) and asymptotic dispersion equation
(30) against a full finite-element treatment. We examine the impact of varying
the aspect ratio $h=2\mathbcal{m}/2\ell=(\bar{b}-\bar{a})/2\bar{\ell}$,
varying the aperture width $\theta_{\mathrm{ap}}$, and varying the filling
fraction $f$ upon the spectral behaviour of the array, as well as its impact
on the Helmholtz resonance frequency. We also evaluate expressions for the
effective inverse density and effective inverse Bulk modulus (compressibility)
for a selection of thin-walled, moderately thick-walled, and specially-scaled
thick-walled Helmholtz resonator arrays in Figures 8 and 9, within the
isotropic (45) and anisotropic (48) approximations given earlier. For
reference we consider wave propagation through air and take $B=0.14183$ MPa
and $\rho=1.2041$ kg/m3 accordingly.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

Figure 3: Comparison of band diagrams (near the $\Gamma$ point) for a two-
dimensional square array of specially-scaled Helmholtz resonators obtained
using multipole methods. Figure LABEL:sub@fig:fold1 shows a spurious result
obtained using the system from Part I (19) with the updated
$\bar{\bar{h}}_{\varepsilon}$ (17), and LABEL:sub@fig:fold2 shows the correct
result obtained using the regularised system (22). Inset: corresponding
fundamental unit cells. Results are given for a dipole truncation $L=1$ and
truncations $L_{\chi_{n}}=L_{F_{m}}=L_{H_{\varepsilon}}=13$ with $\bar{d}=1$,
$\theta_{0}=0$, $\bar{b}=0.3$, $h=100$, and $\theta_{\mathrm{ap}}=\pi/1024$.

In Figure 3 we examine the band diagram for a specially-scaled resonator near
the $\Gamma$ symmetry point, comparing the result obtained from the system in
Part I (with the replacement
$h_{\varepsilon}\mapsto\bar{\bar{h}}_{\varepsilon}$) against that obtained
from the regularised system (22) directly. In the former case, Figure 3a
demonstrates unexpected spectral behaviour in the form of band folding effects
near $k$ values corresponding to $\bar{\bar{h}}_{\varepsilon}\approx 0$. At
the cusps of the folded bands, the group velocity
$\partial\omega/\partial\bar{k}_{Bi}$ is in principle infinite [19], although
as shown in Figure 3b, such features are in fact spurious and are not a
feature of the genuine system. Such behaviour demonstrates the need to
exercise appropriate caution when calculating band diagrams using multipole
methods, although we remark that spurious spectral behaviour can be overcome
using numerical techniques: by searching for both the zero determinant and
vanishing minimum singular value of the matrix system and only considering
those values which satisfy both measures.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

Figure 4: Band diagrams for two-dimensional square arrays of specially-scaled
Helmholtz resonators comparing the finite-element solution (blue curves)
against results obtained using the regularised system (22) for dipole $L=1$
(dashed black curves) and quadrupole $L=3$ (green curves) truncations. Inset:
corresponding fundamental unit cells. Figure LABEL:sub@fig:fold1 corresponds
to $h=100$ ($\bar{a}\approx 0.116$) and $\theta_{\mathrm{ap}}=\pi/1024$, and
LABEL:sub@fig:fold2 corresponds to $\bar{a}=0.1$ and
$\theta_{\mathrm{ap}}=\pi/64$. In both figures, we use the truncations
$L_{\chi_{n}}=L_{F_{m}}=L_{H_{\varepsilon}}=13$ with $\bar{d}=1$,
$\theta_{0}=0$, and $\bar{b}=0.3$.

In Figure 4 we examine the band diagram for two specially-scaled resonator
configurations, comparing results for the regularised system (22) under
different truncations against those obtained using finite element methods. In
general, we find that a dipole approximation (dashed black lines) works quite
well up to the saddle point frequency of the second band surface, with
quadrupole corrections required only for higher frequencies. Hence, in the
figures that follow (Figures 5–7) we consider a quadrupole system truncation
($L=3$) for overall accuracy and a dipole system truncations to derive
asymptotic descriptions. A significant change is observed in the Helmholtz
resonance frequency (equivalently, the maximum frequency of the first band
surface) between the two configurations in Figure 4, which is attributable to
the change in aperture angle. The impact of this parameter is discussed in
further detail below in Figure 10.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

(c) Subfigure 2 list of figures text

(d) Subfigure 3 list of figures text

Figure 5: Band diagrams for a two-dimensional square array of specially-scaled
thick-walled Helmholtz resonators as the channel aspect ratio
$h=\mathbcal{m}/\varepsilon$ is increased: LABEL:sub@fig:thickmultipolesubfig1
$h=10$, LABEL:sub@fig:thickmultipolesubfig2 $h=50$,
LABEL:sub@fig:thickmultipolesubfig3 $h=100$, and
LABEL:sub@fig:thickmultipolesubfig4 $h=150$, with fundamental unit cells
inset. Multipole results from the regularised system (22) are given (green
lines) for a system truncation $L=3$ and truncations
$L_{\chi_{n}}=L_{F_{m}}=L_{H_{\varepsilon}}=13$ with the isotropic
approximation (30) (dashed red lines) superposed, in addition to finite-
element results (blue lines). In the above figures we use $\bar{d}=1$,
$\theta_{0}=0$, $\bar{b}=0.3$, and $\theta_{\mathrm{ap}}=\pi/1024$.

In Figure 5 we compute the band diagrams for a specially-scaled Helmholtz
resonator array as the channel aspect ratio $h$ is varied. Here we observe
that the multipole system, within a quadrupole approximation, is likewise able
to recover the spectral behaviour to excellent accuracy and that the isotropic
approximation (30) is able to recover the first band, determine the width of
the first band gap, and describe the second band surface at low frequencies,
provided that $h$ is not too large. For all values of $h$ we find that the
isotropic approximation works surprisingly well over the range of the first
band surface.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

(c) Subfigure 2 list of figures text

(d) Subfigure 3 list of figures text

Figure 6: Band diagrams for a two-dimensional square array of specially-scaled
thick-walled Helmholtz resonators as the aperture half-angle
$\theta_{\mathrm{a}\mathrm{p}}$ is increased:
LABEL:sub@fig:thicktmultipolesubfig1 $\theta_{\mathrm{a}\mathrm{p}}=\pi/1024$,
LABEL:sub@fig:thicktmultipolesubfig2 $\theta_{\mathrm{a}\mathrm{p}}=\pi/256$,
LABEL:sub@fig:thicktmultipolesubfig3 $\theta_{\mathrm{a}\mathrm{p}}=\pi/64$,
and LABEL:sub@fig:thicktmultipolesubfig4
$\theta_{\mathrm{a}\mathrm{p}}=\pi/16$, with fundamental unit cells inset.
Figure legends and truncation values are identical to those in Figure 5; here
we use $\bar{d}=1$, $\theta_{0}=0$, $\bar{b}=0.3$, and $\bar{a}=0.1$.

Similarly, Figure 6 examines the impact of varying the aperture half-angle
$\theta_{\mathrm{ap}}$ on the band diagram for specially-scaled resonator
arrays. Once more, we find that the multipole system exhibits excellent
performance, and that the isotropic approximation (30) also performs
surprisingly well, provided that the aperture half-angle is not too large. In
fact, for $\theta_{\mathrm{ap}}=\pi/16$ we observe that the first and second
band surfaces are degenerate along $\mathrm{M}\mathrm{Y}$. It is clear that
the Helmholtz resonance/cut-off frequency is considerably sensitive to varying
aperture angle, which we discuss further in Figure 10 below.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

(c) Subfigure 2 list of figures text

(d) Subfigure 3 list of figures text

Figure 7: Band diagrams for a two-dimensional square array of specially-scaled
thick-walled Helmholtz resonators as the filling fraction $f$ is increased:
LABEL:sub@fig:thickfmultipolesubfig1 $f=0.1$,
LABEL:sub@fig:thickfmultipolesubfig2 $f=0.3$,
LABEL:sub@fig:thickfmultipolesubfig3 $f=0.5$, and
LABEL:sub@fig:thickfmultipolesubfig4 $f=0.7$, with fundamental unit cells
inset. Figure legends and truncation values are identical to those in Figure
5; here we use $\bar{d}=1$, $\theta_{0}=0$, $h=100$, and
$\theta_{\mathrm{ap}}=\pi/1024$.

In Figure 7 we determine the band diagrams for a specially-scaled Helmholtz
resonator array, as the filling fraction is varied. As in the preceding
figures, the multipole treatment works remarkably well, even as the outer
resonator wall almost touches the boundaries of the unit cell. This is quite
surprising given that an ever increasing number of multipole orders are
required to compute the band diagram for ideal cylinders in the same limit
[14]. As before, the isotropic approximation (30) is able to recover the first
band surface to exceptional accuracy, and is able to extend into the second
band surface for moderate filling fractions. Even as the filling fraction
approaches the wall-touching limit, the approximation (30) is still able to
determine the width of the first band gap to suitable accuracy.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

(c) Subfigure 2 list of figures text

(d) Subfigure 3 list of figures text

(e) Subfigure 2 list of figures text

(f) Subfigure 3 list of figures text

Figure 8: The effective inverse density and effective inverse Bulk modulus for
a two-dimensional array of thin-walled Helmholtz resonators
LABEL:sub@fig:rhoinvisothin,LABEL:sub@fig:Binvisothin within an isotropic
approximation (45) for $\theta_{0}=0$,
LABEL:sub@fig:rhoinvanisothin,LABEL:sub@fig:Binvanisothin within an
anisotropic approximation (48) for $\theta_{0}=0$, and
LABEL:sub@fig:rhoinvanisothinrot,LABEL:sub@fig:Binvanisothinrot within an
anisotropic approximation (48) for $\theta_{0}=\pi/6$. The shaded red regions
denote all $k$-values above the Helmholtz resonance/cut-off frequency given by
$1-2f_{\varepsilon}\approx 0$ and (49), respectively. Inset: corresponding
fundamental unit cells. In the above figures we use $\bar{d}=1$,
$\theta_{\mathrm{ap}}=\pi/12$, and $\bar{b}=0.3$.

In Figure 8 we examine the effective inverse density and effective inverse
Bulk modulus for an array of thin-walled Helmholtz resonators, specifically
comparing the isotropic and anisotropic expressions (45) and (48). In Figures
8LABEL:sub@fig:rhoinvisothin,LABEL:sub@fig:Binvisothin we observe that within
an isotropic approximation, the band edge coincides with the Helmholtz
resonance frequency ($2\bar{\bar{f}}_{\varepsilon}-1\approx 0$), where the
compressibility likewise diverges. As expected from the representations in
(45), frequency dependence is observed in the inverse Bulk modulus alone. In
Figures 8LABEL:sub@fig:rhoinvanisothin,LABEL:sub@fig:Binvanisothin we observe
an expected frequency dependence in one component of the effective inverse
density matrix as well as for the inverse Bulk modulus, when anisotropy is
considered, see (48). Here,
$(\rho_{\mathrm{eff}})^{-1}_{\mathrm{y}\mathrm{y}}$ also diverges at the
Helmholtz resonance frequency (49). Similarly, in Figures
8LABEL:sub@fig:rhoinvanisothinrot,LABEL:sub@fig:Binvanisothinrot, we observe
frequency dependence and divergence at the resonance frequency (49) for all
entries in $(\rho_{\mathrm{eff}})^{-1}_{ij}$ as well as for the inverse Bulk
modulus, when $\theta_{0}=\pi/6$. By rotating the resonator we have reduced
the symmetry properties of the medium, allowing for stronger anisotropy and
dispersion. In any case, despite the emergence of anisotropy in the effective
inverse density, no changes in sign are observed in the coefficients over the
range of the first band surface (
$(\rho_{\mathrm{eff}})^{-1}_{\mathrm{x}\mathrm{y}}$ is either zero or
negative). This suggests that exotic effects such as negative refraction are
not supported on the first band surface for Helmholtz resonator arrays,
however we expect that negative refraction is supported at higher frequencies
where the correct band curvature is exhibited [20].

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

(c) Subfigure 2 list of figures text

(d) Subfigure 3 list of figures text

Figure 9: The effective inverse density and effective inverse Bulk modulus
LABEL:sub@fig:rhoinvisothin2,LABEL:sub@fig:Binvisothin2 for a two-dimensional
array of moderately thick-walled Helmholtz resonators within an anisotropic
approximation (48);
LABEL:sub@fig:rhoinvanisothin2,LABEL:sub@fig:Binvanisothin2 for an array of
specially-scaled Helmholtz resonators within an isotropic approximation (45).
Inset: corresponding fundamental unit cells. In
LABEL:sub@fig:rhoinvisothin2,LABEL:sub@fig:Binvisothin2 the shaded red region
denotes all $k$-values above the Helmholtz resonance/cut-off frequency (49),
where we use $\bar{d}=1$, $\theta_{0}=0$, $\theta_{\mathrm{ap}}=\pi/48$,
$h=3$, and $\bar{b}=0.3$. In
LABEL:sub@fig:rhoinvanisothin2,LABEL:sub@fig:Binvanisothin2 the shaded red
region denotes the width of the first band gap as shown in Figure
5LABEL:sub@fig:thickmultipolesubfig3, where we use $\bar{d}=1$,
$\theta_{0}=0$, $\theta_{\mathrm{ap}}=\pi/1024$, $h=100$, and $\bar{b}=0.3$.

In Figures 9LABEL:sub@fig:rhoinvisothin2,LABEL:sub@fig:Binvisothin2 we
consider the effective inverse density and inverse Bulk modulus functions,
within an anisotropic description, corresponding to a moderately thick-walled
resonator array. These curves are obtained via (48) with the replacement
$f_{\varepsilon}\mapsto\check{f}_{\varepsilon}$ (results for the isotropic
descriptions are not included here as these are qualitatively similar to those
given for the thin-walled case). As observed in the thin-walled case, we see
frequency dependence in both material tensors over the frequency range of the
first band surface, which is accompanied by unexpected singular behaviour
above the Helmholtz resonance frequency. No firm conclusions may be drawn from
this, however, as this behaviour lies outside the region of validity for the
expressions. In Figures
9LABEL:sub@fig:rhoinvanisothin2,LABEL:sub@fig:Binvanisothin2 we consider a
specially-scaled resonator array (i.e., (45) with the replacement
$f_{\varepsilon}\mapsto\bar{\bar{f}}_{\varepsilon}$) which exhibits much more
interesting behaviour. In particular, we observe a pole in
$B_{\mathrm{eff}}^{-1}$ at low frequencies, where the width of the first band
gap may be defined as the interval between the Helmholtz resonance frequency
given by $1-2\bar{\bar{f}}_{\varepsilon}\approx 0$ and the zero of
$B_{\mathrm{eff}}^{-1}$; that is, the first band gap corresponds to the
interval where the effective bulk modulus is negative. Once this physical
parameter returns to positive values we find that these isotropic descriptions
extend well into the range of the second band surface, as seen in Figure
5LABEL:sub@fig:thickmultipolesubfig3. For reference, we estimate our
descriptions (45) to hold over the approximate range $0\leq k\leq 1.5$ for
this example, which supports the assertion that the isotropic description for
specially-scaled resonator arrays appears to be valid over a significantly
broader frequency range than for the resonators discussed in Part I.

Finally, in Figure 10, to complement the band diagram figures outlined in
Figures 5 and 6, we solve $2\bar{\bar{f}}_{\varepsilon}-1\approx 0$ to track
the Helmholtz resonance/cut-off frequency as the aspect ratio $h$ and aperture
half-angle $\theta_{\mathrm{ap}}$ are varied. In Figure 10a we observe that
the cut-off frequency for the first band surface is able to achieve a minimum
of $k\approx 0.369$ at $h\approx 52$ for a fixed
$\theta_{\mathrm{ap}}=\pi/1024$, and that as the aspect ratio $h$ becomes
large, the cut-off frequency grows larger with polynomial scaling. The
possibility of achieving a cutoff frequency minimum, by tuning the wall
thickness, may prove useful for those involved in the design of acoustic
metamaterials. Likewise, Figure 10b demonstrates that the resonant frequency
decreases with decreasing aperture half-angle, as expected, and scales as
$k_{\mathrm{max}}\propto\theta_{\mathrm{ap}}^{1/2}$ as
$\theta_{\mathrm{ap}}\rightarrow 0$.

(a) Subfigure 6 list of figures text

(b) Subfigure 1 list of figures text

Figure 10: Tracking the cut-off/Helmholtz resonance frequency (by solving
$2\bar{\bar{f}}_{\varepsilon}-1\approx 0$) in the specially-scaled setting as
LABEL:sub@fig:helmtrack1 the aspect ratio of the resonator neck
$h={\mathbcal{m}}/{\ell}$ is varied, as shown in Figure 5, and as
LABEL:sub@fig:helmtrack2 the aperture half-angle $\theta_{\mathrm{ap}}$ is
varied, as demonstrated in Figure 6.

##  8 Discussion

In this paper, we have presented a matched asymptotic-multipole procedure for
determining the band structure of an acoustic metamaterial comprising a two-
dimensional array of Helmholtz resonators that possess large wall thicknesses
and narrow neck widths. We have also derived a compact dispersion equation
which is able to describe the first band surface, first band gap, and
frequencies well into the second band surface, over a range of practical
settings. In addition, we have outlined a homogenisation procedure for thin-
walled, moderately thick-walled, and specially-scaled resonator arrays,
presenting analytical forms for the effective inverse density tensor and the
effective inverse Bulk modulus. The effective response functions derived
extend well beyond the quasistatic limit to higher frequencies, depending on
the resonator geometry. Furthermore, we have also derived closed-form
representations for the Helmholtz resonance frequency in the specially-scaled
setting.

We demonstrate that the specially-scaled resonators are able to achieve
extremely low Helmholtz resonance frequencies, in contrast to thin- and
moderately-thick walled resonators presented in Part I, which has a marked
impact on the performance of the array. By incorporating long neck widths
$2\bar{\mathbcal{m}}$ into the formulation we provide an additional degree of
freedom for controlling the frequency range of the first band surface,
indirectly controlling features such as the low-frequency phase and group
velocity. Incidentally, we remark that the Helmholtz resonance frequency is
often approximated in the literature via the form [21, Eq. (5.3.12)]

$\omega_{\mathrm{max}}\approx\sqrt{\frac{B}{\rho}}\sqrt{\frac{A}{LV}},\quad\mbox{
or equivalently},\quad
k_{\mathrm{max}}\approx\frac{1}{\bar{a}}\sqrt{\dfrac{1}{\pi\bar{\mathbcal{m}}/\bar{\ell}}}$
(50)

where $A$ is the total aperture width $2\bar{\ell}$, $L$ denotes the length of
the resonator neck $2\bar{\mathbcal{m}}$, and $V$ denotes the enclosed
resonator area $\pi\bar{a}^{2}$. The major disadvantage of (50) is that the
approximation is quite crude, as it treats the neck as distinct from the
enclosed volume, and as a result frequently requires correction factors and an
effective neck length to recover accuracy. We stress that our Helmholtz
resonance expressions derived in (35b), (38), and in Part I, do not require
any such corrections, and apply to both single resonators and arrays of
resonators alike, based purely on asymptotic scale separation via the small
parameter $\varepsilon$.

On comparing results from Parts I and II, we find that anisotropy in both the
band structure and in the effective tensors is greatest for thin-walled
resonators, with anisotropic effects considerably reduced as we approach the
specially-scaled limit. In fact, we find that specially-scaled resonators may
be treated as almost-isotropic media. Note however that results from Part I
are not recovered as $\bar{\mathbcal{m}}\rightarrow 0$ within the current
formulation, or vice versa, due to the different underlying assumptions of the
inner solutions; hence the reason for treating them separately. A core
advantage of our matched asymptotics-multipole treatment is that it avoids the
need for extensive fully-numerical procedures, such a finite-element methods,
for determining the low-frequency band structure of acoustic metamaterials.
Fully-numerical procedures require intensive meshing/sampling inside the neck
region as it becomes increasingly thin, and although the computational domain
is two-dimensional, meshing requirements can massively increase computation
times and resource requirements. In contrast, our asymptotic dispersion
equations provide rapidly-evaluable closed-form representations for band
surfaces over a wide frequency range. It is of interest to apply the
techniques outlined here to related geometries, such as resonators with long
maze-like channels [22] or resonators nested within resonators [23, 24], and
this is presently under investigation by the authors.

In relation to our homogenisation treatment, we emphasise that composite
materials and metamaterials have the potential to exhibit generalised
constitutive relations [25, 26]. Such behaviour occurs widely across the
acoustics and elasticity literature, where these materials are known as Willis
media, as well as in the electromagnetics literature, where they are known as
bi-isotropic or bi-anisotropic media [27]. A significant body of literature
has developed in recent years with such effects in mind [28, 29, 30, 31]. It
would be of interest to examine Willis coupling effects in resonator arrays in
further detail. Recently, a body of experimental and theoretical work has
emerged demonstrating analogues to Willis coupling in wave scattering by a
single thick-walled Helmholtz resonator [22, 32] which could offer interesting
potential applications.

Finally, there are points regarding passivity and causality in resonator
arrays, and the behaviour of matched asymptotic expansion solutions, which are
worth investigating and this will be reported on by the authors in a
forthcoming article.

##  Appendix A Asymptotic forms of the lattice sums

As discussed in the Appendix of Part I, the lattice sums $S_{m}^{\mathrm{Y}}$
emerge frequently in the study of periodic media, but are conditionally
convergent in their most direct form and so must be regularised in order to
retrieve physically meaningful results. The convergent expressions given in
Part I are used in all relevant numerical computations, whereas for the
purposes of asymptotic analysis we follow the procedures outlined in [33, 34],
to derive and present leading terms for the first few lattice sums of a square
array in the form:

$\begin{aligned} S_{0}^{\mathrm{Y}}&\sim-\frac{4f}{\pi
b^{2}\left(k_{\mathrm{B}}^{2}-1\right)}+\frac{1}{\pi}\left[-2\gamma_{\mathrm{e}}+\log\left(\frac{f\Gamma(\tfrac{1}{4})^{4}}{\pi^{2}b^{2}}\right)\right],\\\
S_{1}^{\mathrm{Y}}&\sim-\frac{4\mathrm{i}fk_{\mathrm{B}}}{\pi
b^{2}\left(k_{\mathrm{B}}^{2}-1\right)}\mathrm{e}^{\mathrm{i}\theta_{\mathrm{B}}}+\left(\frac{\mathrm{i}k_{\mathrm{B}}}{\pi}\right)\mathrm{e}^{\mathrm{i}\theta_{\mathrm{B}}},\\\
S_{2}^{\mathrm{Y}}&\sim\frac{4fk_{\mathrm{B}}^{2}}{\pi
b^{2}\left(k_{\mathrm{B}}^{2}-1\right)}\mathrm{e}^{2\mathrm{i}\theta_{\mathrm{B}}}+\left[\frac{k_{\mathrm{B}}^{2}\Gamma(\frac{1}{4})^{8}}{384\pi^{5}}\mathrm{e}^{-2\mathrm{i}\theta_{\mathrm{B}}}-\frac{k_{\mathrm{B}}^{2}}{2\pi}\mathrm{e}^{2\mathrm{i}\theta_{\mathrm{B}}}\right],\\\
S_{3}^{\mathrm{Y}}&\sim\frac{1}{b^{2}}\left[\frac{4\mathrm{i}fk_{\mathrm{B}}^{3}}{\pi\left(k_{\mathrm{B}}^{2}-1\right)}\mathrm{e}^{3\mathrm{i}\theta_{\mathrm{B}}}-\frac{\mathrm{i}k_{\mathrm{B}}f\Gamma(\tfrac{1}{4})^{8}}{48\pi^{5}}\mathrm{e}^{-\mathrm{i}\theta_{\mathrm{B}}}\right],\\\
S_{4}^{\mathrm{Y}}&\sim-\frac{f^{2}\Gamma(\frac{1}{4})^{8}}{10\pi^{5}b^{4}},\end{aligned}\qquad\begin{aligned}
S_{5}^{\mathrm{Y}}&\sim-\frac{2\mathrm{i}f^{2}k_{\mathrm{B}}\Gamma(\tfrac{1}{4})^{8}}{5\pi^{5}b^{4}}\mathrm{e}^{\mathrm{i}\theta_{\mathrm{B}}},\\\
S_{6}^{\mathrm{Y}}\sim\frac{1}{b^{4}}\left[\frac{f^{2}k_{\mathrm{B}}^{2}\Gamma(\tfrac{1}{4})^{8}}{\pi^{5}}\mathrm{e}^{2\mathrm{i}\theta_{\mathrm{B}}}\right.&+\left.\frac{5f^{2}k_{\mathrm{B}}^{2}\Gamma(\tfrac{1}{4})^{16}}{960\pi^{9}}\mathrm{e}^{-2\mathrm{i}\theta_{\mathrm{B}}}\right],\\\
S_{7}^{\mathrm{Y}}&\sim-\frac{\mathrm{i}f^{3}k_{\mathrm{B}}\Gamma(\tfrac{1}{4})^{16}}{14\pi^{9}b^{6}}\mathrm{e}^{-\mathrm{i}\theta_{\mathrm{B}}},\\\
S_{8}^{\mathrm{Y}}&\sim-\frac{3f^{4}\Gamma(\tfrac{1}{4})^{16}}{5\pi^{9}b^{8}},\\\
S_{9}^{\mathrm{Y}}&\sim-\frac{24\mathrm{i}\mathrm{e}^{\mathrm{i}\theta_{\mathrm{B}}}f^{4}k_{\mathrm{B}}\Gamma(\tfrac{1}{4})^{16}}{5\pi^{9}b^{8}}.\end{aligned}$
(51)

For reference, we remark that these sums are strictly defined in terms of the
lattice period $d$. However, we have introduced $f=\pi b^{2}/d^{2}$, the
cylinder area fraction, here since the limit of $b\rightarrow 0$ (with $f$
fixed) is considered in our analysis.

## Acknowledgements

I.D.A. acknowledges support from a Royal Society Industry Fellowship. This
work was also supported by EPSRC grant no EP/R014604/1 whilst I.D.A. held the
position of Director of the Isaac Newton Institute Cambridge.

## References

  * [1] Schweizer B. 2017 Resonance meets homogenization. J. Dtsch. Math.-Verein. 119, 31–51.
  * [2] Lamacz A, Schweizer B. 2016 Effective acoustic properties of a meta-material consisting of small Helmholtz resonators. arXiv preprint arXiv:1603.05395.
  * [3] Kohn RV, Shipman SP. 2008 Magnetism and homogenization of microresonators. Multiscale Model. Sim. 7, 62–92.
  * [4] Bouchitté G, Felbacq D. 2004 Homogenization near resonances and artificial magnetism from dielectrics. C. R. Math. 339, 377–382.
  * [5] Zhikov V. 2005 On spectrum gaps of some divergent elliptic operators with periodic coefficients. St. Petersb. Math. J. 16, 773–790.
  * [6] Haberman MR, Guild MD. 2016 Acoustic metamaterials. Phys. Today 69, 42–48.
  * [7] Howe M. 1976 On the Helmholtz resonator. J. Sound. Vib. 45, 427–440.
  * [8] Von Helmholtz H. 1860 Theorie der Luftschwingungen in Röhren mit offenen Enden. Crelle pp. 1–72.
  * [9] Rayleigh L. 1870 On the theory of resonance. Phil. Trans. Roy. Soc. Lond. pp. 77–118.
  * [10] Schelkunoff SA, Friis HT. 1952 Antennas: theory and practice. NY: John Wiley and Sons.
  * [11] Crighton DG, Dowling AP, Ffowcs-Williams JE, Heckl M, Leppington FG. 1992 Modern methods in analytical acoustics lecture notes. Berlin: Springer-Verlag.
  * [12] Cotterill PA, Parnell WJ, Abrahams ID, Miller R, Thorpe M. 2015 The time-harmonic antiplane elastic response of a constrained layer. J. Sound Vib. 348, 167–184.
  * [13] Sloane NJA, Inc. TOF. 2021 The on-line encyclopedia of integer sequences. .
  * [14] Movchan AB, Movchan NV, Poulton CG. 2002 Asymptotic models of fields in dilute and densely packed composites. London, UK: Imperial College Press.
  * [15] Torrent D, Sánchez-Dehesa J. 2006 Effective parameters of clusters of cylinders embedded in a nonviscous fluid or gas. Phys. Rev. B 74, 224305.
  * [16] Martin P, Maurel A, Parnell W. 2010 Estimating the dynamic effective mass density of random composites. J. Acoust. Soc. Am. 128, 571–577.
  * [17] Llewellyn-Smith SG, Davis AMJ. 2010 The split ring resonator. Proc. Roy. Soc. A 466, 3117–3134.
  * [18] Norris AN. 2015 Acoustic cloaking. Acoust. Today 11, 38–46.
  * [19] Chen PY, Poulton C, Asatryan A, Steel M, Botten LC, De Sterke CM, McPhedran R. 2011 Folded bands in metamaterial photonic crystals. N. J. Phys. 13, 053007.
  * [20] Smith MJA, McPhedran RC, Poulton CG, Meylan MH. 2012 Negative refraction and dispersion phenomena in platonic clusters. Waves Rand. Comp. Media 22, 435–458.
  * [21] Howe MS. 1998 Acoustics of fluid-structure interactions. Cambridge, UK: Cambridge University Press.
  * [22] Quan L, Ra’di Y, Sounas DL, Alù A. 2018 Maximum Willis coupling in acoustic scatterers. Phys. Rev. Lett. 120, 254301.
  * [23] Elford DP, Chalmers L, Kusmartsev FV, Swallowe GM. 2011 Matryoshka locally resonant sonic crystal. J. Acoust. Soc. Am. 130, 2746–2755.
  * [24] Montiel F, Chung H. 2020 Planar acoustic scattering by a multi-layered split ring resonator. J. Acoust. Soc. Am. 148, 3698–3708.
  * [25] Milton GW, Willis JR. 2007 On modifications of Newton’s second law and linear continuum elastodynamics. Proc. Roy. Soc. A 463, 855–880.
  * [26] Willis JR. 2011 Effective constitutive relations for waves in composites and metamaterials. Proc. Roy. Soc. A 467, 1865–1879.
  * [27] Lindell IV, Sihvola A, Viitanen A, Tretyakov S. 1994 Electromagnetic waves in chiral and bi-isotropic media. Norwood MA, USA: Artech House.
  * [28] Norris AN, Shuvalov A, Kutsenko A. 2012 Analytical formulation of three-dimensional dynamic homogenization for periodic elastic systems. Proc. Roy. Soc. A 468, 1629–1651.
  * [29] Torrent D, Pennec Y, Djafari-Rouhani B. 2015 Resonant and nonlocal properties of phononic metasolids. Phys. Rev. B 92, 174110.
  * [30] Muhlestein MB, Sieck CF, Alù A, Haberman MR. 2016 Reciprocity, passivity and causality in Willis materials. Proc. Roy. Soc. A 472, 20160604\.
  * [31] Muhlestein MB, Sieck CF, Wilson PS, Haberman MR. 2017 Experimental evidence of Willis coupling in a one-dimensional effective material element. Nat. Comm. 8, 1–9.
  * [32] Melnikov A, Chiang YK, Quan L, Oberst S, Alù A, Marburg S, Powell D. 2019 Acoustic meta-atom with experimentally verified maximum Willis coupling. Nat. Comm. 10, 1–7.
  * [33] McPhedran R, Poulton C, Nicorovici N, Movchan A. 1996 Low frequency corrections to the static effective dielectric constant of a two-dimensional composite material. Proc. Roy. Soc. A 452, 2231–2245.
  * [34] Chen PY, Smith MJA, McPhedran RC. 2018 Evaluation and regularization of phase-modulated Eisenstein series and application to double Schlömilch-type sums. J. Math. Phys. 59, 072902.


* [60] E. Ackerman, “Meet the Lunar Gateway’s Robot Caretakers: With people seldom on board, the new space station will rely on autonomy,” _IEEE Spectrum_ , Apr. 2022. [Online]. Available: https://spectrum.ieee.org/lunar-gateway-robots
  * [61] S. Paul and K. Huper, “An analog circuit for eigenvalue calculation and rank filtering,” _IEEE Trans. Circuits Syst. I. Fundam. Theory Appl._ , vol. 41, no. 11, pp. 736–740, Nov. 1994.
  * [62] S. Suh, A. Basu, C. Schlottmann, P. E. Hasler, and J. R. Barry, “Low-power discrete Fourier transform for OFDM: A programmable analog approach,” _IEEE Trans. Circuits Syst. I, Reg. Papers_ , vol. 58, no. 2, pp. 290–298, Feb. 2011.
  * [63] M. Gu and S. Chakrabartty, “Synthesis of bias-scalable CMOS analog computational circuits using margin propagation,” _IEEE Trans. Circuits Syst. I, Reg. Papers_ , vol. 59, no. 2, pp. 243–254, Feb. 2012.
  * [64] W. Haensch, T. Gokmen, and R. Puri, “The next generation of deep learning hardware: Analog computing,” _Proc. IEEE_ , vol. 107, no. 1, pp. 108–122, Jan. 2019.
  * [65] N. Udayanga, S. I. Hariharan, S. Mandal, L. Belostotski, L. T. Bruton, and A. Madanayake, “Continuous-time algorithms for solving Maxwell’s equations using analog circuits,” _IEEE Trans. Circuits Syst. I, Reg. Papers_ , vol. 66, no. 10, pp. 3941–3954, Oct. 2019.
  * [66] O. Bournez and A. Pouly, “A survey on analog models of computation,” in _Handbook of Computability and Complexity in Analysis_ , ser. Theory and Applications of Computability, V. Brattka and P. Hertling, Eds. Springer, 2021, pp. 173–226.
  * [67] A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” _Science_ , vol. 343, no. 6167, pp. 160–163, 2014.
  * [68] F. Zangeneh-Nejad and R. Fleury, “Topological analog signal processing,” _Nat. Commun._ , vol. 10, no. 2058, pp. 1223–1231, May 2019.
  * [69] C. E. Shannon, “Mathematical theory of the differential analyzer,” _Journal of Mathematics and Physics_ , vol. 20, no. 1–4, pp. 337–354, Apr. 1941.
  * [70] D. S. Graça and J. F. Costa, “Analog computers and recursive functions over the reals,” _J. Complexity_ , vol. 19, no. 5, pp. 644–664, Oct. 2003\.
  * [71] V. Brattka, P. Hertling, and K. Weihrauch, “A Tutorial on Computable Analysis,” in _New Computational Paradigms_ , S. B. Cooper, B. Löwe, and A. Sorbi, Eds. New York, USA: Springer Science+Business Media, LLC, 2008.
  * [72] W. Rudin, _Real and Complex Analysis_ , 3rd ed. Boston, etc.: McGraw-Hill, 1987.