Topologically protected edge modes in one-dimensional chains of subwavelength resonators

TL;DR

This paper develops a design for subwavelength resonator chains that supports topologically protected edge modes, enabling robust wave localization at deep subwavelength scales despite imperfections.

Contribution

It introduces a topological framework for designing stable, subwavelength waveguiding structures using high-contrast resonator chains.

Findings

Topologically protected edge modes are demonstrated in finite resonator chains.

Wave localization is robust against random imperfections.

Design principles for stable subwavelength waveguides are established.

Abstract

The goal of this paper is to advance the development of wave-guiding subwavelength crystals by developing designs whose properties are stable with respect to imperfections in their construction. In particular, we make use of a locally resonant subwavelength structure, composed of a chain of high-contrast resonators, to trap waves at deep subwavelength scales. We first study an infinite chain of subwavelength resonator dimers and define topological quantities that capture the structure's wave transmission properties. Using this for guidance, we design a finite crystal that is shown to have wave localization properties, at subwavelength scales, that are robust with respect to random imperfections.

Equations

$$G^k(x,y):=-\frac{e^{\mathrm{i}k|x-y|}}{4\pi |x-y|},x,y\in {\mathbb{R}}^3,k\ge 0.$$

$${\mathcal{S}}_D^k[\varphi ](x):={\int }_{\partial D}{G}^k(x,y)\varphi (y)\mathrm{d}\sigma (y),x\in {\mathbb{R}}^3.$$

$${\mathcal{K}}_D^{k,\ast }[\varphi ](x):={\int }_{\partial D}\frac{\partial }{\partial {\nu }_x}{G}^k(x,y)\varphi (y)\mathrm{d}\sigma (y),x\in \partial D,$$

$$\mathcal{S}_{D}^{k}[\phi]\big{|}_{+}=\mathcal{S}_{D}^{k}[\phi]\big{|}_{-},$$

$$\frac{\partial}{\partial\nu}\mathcal{S}_{D}^{k}[\phi]\Big{|}_{\pm}=\left(\pm\frac{1}{2}I+\mathcal{K}_{D}^{k,*}\right)[\phi],$$

$$\mathcal{F}[f](x,\alpha ):=\sum _{m\in \mathbb{Z}}f(x-Lm)e^{\mathrm{i}L\alpha m}.$$

$${\mathcal{F}}^{-1}[g](x)=\frac{1}{2\pi }{\int }_{Y^{\ast }}g(x,\alpha )\mathrm{d}\alpha ,x\in \mathbb{R},$$

$$G^{\alpha ,k}(x,y):=-\sum _{m\in \mathbb{Z}}\frac{e^{\mathrm{i}k|x-y-(Lm,0,0)|}}{4\pi |x-y-(Lm,0,0)|}{e}^{\mathrm{i}\alpha Lm}.$$

$${\mathcal{S}}_D^{\alpha ,k}[\varphi ](x):={\int }_{\partial D}{G}^{\alpha ,k}(x,y)\varphi (y)\mathrm{d}\sigma (y),x\in {\mathbb{R}}^3.$$

$$\mathcal{S}_{D}^{\alpha,k}[\phi]\big{|}_{+}=\mathcal{S}_{D}^{\alpha,k}[\phi]\big{|}_{-},$$

$$\frac{\partial}{\partial\nu}\mathcal{S}_{D}^{\alpha,k}[\phi]\Big{|}_{\pm}=\left(\pm\frac{1}{2}I+(\mathcal{K}_{D}^{-\alpha,k})^{*}\right)[\phi]\quad\mbox{on}\leavevmode\nobreak\ \partial D,$$

$$({\mathcal{K}}_D^{-\alpha ,k}{)}^{\ast }[\varphi ](x):={\int }_{\partial D}\frac{\partial }{\partial {\nu }_x}{G}^{\alpha ,k}(x,y)\varphi (y)\mathrm{d}\sigma (y).$$

$$\sigma (\mathcal{L})=\bigcup _{\alpha \in Y^{\ast }}\sigma (\mathcal{L}(\alpha )).$$

$$A_n:=\mathrm{i}⟨u_n^{\alpha },\frac{\partial }{\partial \alpha }{u}_n^{\alpha }⟩.$$

$$\theta ={\int }_{{\alpha }_1}^{{\alpha }_2}{A}_n\mathrm{d}\alpha .$$

$$\varphi_{n}^{z}:=\mathrm{i}\mkern 1.0mu\int_{Y^{*}}\big{\langle}u_{n}^{\alpha},\frac{\partial}{\partial\alpha}u_{n}^{\alpha}\big{\rangle}\>\mathrm{d}\alpha,$$

$$R_1{D}_1=D_1,R_2{D}_2=D_2,$$

$$D=-D.$$

$$\mathcal{C}:=\bigcup _{m\in \mathbb{Z}}\left(D+(mL,0,0)\right).$$

$$\left\{\begin{array}[]{ll}\displaystyle\Delta{u}+\omega^{2}{u}=0&\text{in }\mathbb{R}^{3}\setminus\partial\mathcal{C},\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle{u}|_{+}-{u}|_{-}=0&\text{on }\partial\mathcal{C},\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\delta\frac{\partial{u}}{\partial\nu}\bigg{|}_{+}-\frac{\partial{u}}{\partial\nu}\bigg{|}_{-}=0&\text{on }\partial\mathcal{C},\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle u(x_{1},x_{2},x_{3})&\text{satisfies the outgoing radiation condition as }\sqrt{x_{2}^{2}+x_{3}^{2}}\rightarrow\infty.\end{array}\right.$$

$$\delta \ll 1.$$

$$\mathcal{L}:=-\frac{1}{\delta+(1-{\delta})\chi_{\mathcal{C}}}\nabla\cdot\left(\Big{(}\delta+(1-{\delta})\chi_{\mathcal{C}}\Big{)}\nabla\right),$$

$$\left\{\begin{array}[]{ll}\displaystyle\Delta u_{\alpha}+\omega^{2}{u_{\alpha}}=0&\text{in }\mathbb{R}^{3}\setminus\partial\mathcal{C},\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle{u_{\alpha}}|_{+}-{u_{\alpha}}|_{-}=0&\text{on }\partial\mathcal{C},\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle\delta\frac{\partial{u_{\alpha}}}{\partial\nu}\bigg{|}_{+}-\frac{\partial{u_{\alpha}}}{\partial\nu}\bigg{|}_{-}=0&\text{on }\partial\mathcal{C},\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle e^{-\mathrm{i}\mkern 1.0mu\alpha_{1}x_{1}}u_{\alpha}(x_{1},x_{2},x_{3})&\mbox{is periodic in }x_{1},\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle u_{\alpha}(x_{1},x_{2},x_{3})&\text{satisfies the $\alpha$-quasiperiodic outgoing radiation condition}\\ &\hskip 14.22636pt\text{as }\sqrt{x_{2}^{2}+x_{3}^{2}}\rightarrow\infty.\end{array}\right.$$

$$\mathcal{L}(\alpha):=-\frac{1}{\delta+(1-{\delta})\chi_{\mathcal{C}}}\nabla\cdot\left(\Big{(}\delta+(1-{\delta})\chi_{\mathcal{C}}\Big{)}\nabla\right),$$

$$\sigma =\bigcup _{\alpha \in Y^{\ast }}{\sigma }_{\alpha },$$

$$u_{\alpha }={\mathcal{S}}_D^{\alpha ,\omega }[{\Psi }^{\alpha }],$$

$${\mathcal{A}}^{\alpha }(\omega ,\delta )[{\Psi }^{\alpha }]=0,$$

$${\mathcal{A}}^{\alpha }(\omega ,\delta ):=-\lambda I+{\left({\mathcal{K}}_D^{-\alpha ,\omega }\right)}^{\ast },\lambda :=\frac{1+\delta }{2(1-\delta )}.$$

$$\{\begin{array}{ll}\Delta V_j^{\alpha }=0& \text{mbox}inY\setminus D,\\ V_j^{\alpha }={\delta }_{ij}& \text{mbox}on\partial D_i,\\ V_j^{\alpha }(x+(mL,0,0))=e^{\mathrm{i}\alpha m}{V}_j^{\alpha }(x)& \forall m\in \mathbb{Z},\\ V_j^{\alpha }(x_1,x_2,x_3)=O\left(\frac{1}{\sqrt{x_2^2+x_3^2}}\right)& \text{as }\sqrt{x_2^2+x_2^2}\to \infty ,\text{ uniformly in }{x}_1,\end{array}$$

$$C_{ij}^{\alpha }:={\int }_{Y\setminus D}\overline{\nabla V_i^{\alpha }}\cdot \nabla V_j^{\alpha }\mathrm{d}x,i,j=1,2.$$

Full text

Mathematics Subject Classification (MSC2000): 35J05, 35C20, 35P20.

Keywords: subwavelength resonance, subwavelength phononic and photonic crystals, topological nanomaterials, edge states.

1 Introduction

The ability to manipulate and guide the propagation of waves on subwavelength scales is important for many different physical applications. In the contexts of nanophotonics and nanophononics, subwavelength crystalline structures have, in particular, been shown to have desirable properties. Here, subwavelength means that the size of the repeating unit cell is several magnitudes smaller than the incident wavelengths. It was recently shown, for example, that one can design subwavelength crystals with low ranges of frequencies that cannot propagate (known as subwavelength band gaps) [6] and can localize (or trap) specific frequencies at subwavelength scales by introducing a defect [3]. However, one limitation of such designs is that their properties are often very sensitive to imperfections in the crystal’s structure. In order to be able to feasibly manufacture wave-guiding devices, it is important that we are able to design subwavelength crystals that exhibit stability with respect to geometric errors.

We take inspiration from quantum mechanics. So-called topological insulators have been extensively studied for their electronic properties, in the setting of the Schrödinger operator [15, 17, 20, 21, 22, 23, 32]. The principle that underpins the design of these structures is that one is able to define topological invariants which capture the crystal’s wave propagation properties. Then, if part of a crystalline structure is replaced with an arrangement that is associated with a different value of this invariant, not only will certain frequencies be localized to the interface (as predicted by the classical theory for crystals with defects) but this behaviour will be stable with respect to imperfections. These eigenmodes are known as edge modes and we say that they are topologically protected to refer to their robustness.

One of the most classical examples from quantum mechanics is the well-known Su-Schrieffer-Heeger (SSH) model [42]. Originally introduced to study the electrical properties of polyacetylene, the SSH model consists of a chain of atoms arranged as dimers (similar to that depicted in Figure 1). In the case of one-dimensional crystals such as this, the natural choice of topological invariant is the Zak phase [47]. Qualitatively, a non-zero Zak phase means that the crystal has undergone band inversion, meaning that at some point in the Brillouin zone the monopole/dipole nature of the first/second Bloch eigenmodes has swapped. In this way, the Zak phase captures the crystal’s wave propagation properties. The Zak phased was measured experimentally by [13], and in the SSH model this can take two discrete values depending on whether the atoms in each dimer are closer to each other than they are to the next dimer in the chain. In higher dimensional crystals, topological indices are similarly dependent on the symmetry of the crystals [40]. If one takes two SSH chains with different Zak phases and joins half of one chain to half of the other to form a new crystal, this crystal will exhibit a topologically protected edge mode at the interface. This principle is known as the bulk-boundary correspondence in quantum settings [15, 16, 26, 27, 28, 35, 39]. Here, the term bulk is used to refer to parts of a crystal that are away from an edge (and so resemble an infinite, defect-free crystalline structure).

Understanding why topologically protected edge modes are stable to local perturbations is subtle, and doing so precisely is very much an open question. It can be argued that, due to (chiral) symmetry, these crystals not only have band gaps but the frequencies associated with edge modes occur in the centre of the band gap. We call them midgap frequencies. This is in sharp contrast to conventional, unprotected defect frequencies, which typically emerge from the edge of the band gap [3]. As a result, a small imperfection in the subwavelength crystal will not be able to move a topologically protected frequency out of the band gap, while an unprotected frequency is often lost amongst the bulk frequencies. Moreover, if the perturbation preserves the crystal’s symmetry, the frequency of the edge mode will be very stable, experiencing much smaller variations compared to the other subwavelength resonant frequencies. These effects are typically explained as a consequence of the different topological properties on either side of the edge (see, for example, [36] for a review of topological phases in acoustic systems).

Subwavelength topological photonic and phononic crystals, based on locally resonant crystalline structures with large material contrasts, have been studied both numerically and experimentally in [41, 43, 44, 45, 46]. Subwavelength crystals allow for the manipulation and localization of waves on very small spatial scales and are therefore very useful in physical applications, especially situations where the operating wavelengths are very large. Recently, topological properties of acoustic waves in SSH chains and honeycomb lattices of subwavelength resonators have been numerically and experimentally explored [11, 48, 50, 51]. In this work, we study a subwavelength crystal exhibiting a topologically non-trivial band gap. The crystal consists of a chain of subwavelength resonators arranged as dimers, similar to the SSH model (see Figures 1, LABEL: and 6). Wave propagation in the resonant structure is modelled by a high-contrast Helmholtz problem. High material contrasts are an essential prerequisite for the existence of resonant behaviours on subwavelength scales [2, 37]. Such problems arise naturally in the context of both nanophotonic and nanophononic structures [1, 2, 4]. Around this frequency, a single resonator in free-space scatters waves with a greatly enhanced amplitude. If a second resonator is introduced, coupling interactions will occur giving a system that has both monopole and dipole resonance modes [7]. This pattern continues for larger systems [1].

Initially, in Section 3, we set out to study the bulk properties of an infinitely periodic chain of subwavelength resonator dimers. Using Floquet-Bloch theory, we are able to analytically derive the resonant frequencies and associated eigenmodes of this crystal, and further prove that there exists a non-trivial band gap. Motivated by the use of the Zak phase in quantum mechanics, as well as the work of [33, 36, 38, 49] in photonics and phononics, we define a topological invariant which we will also refer to as the Zak phase. We prove that the Zak phase takes different values for different geometries and in the dilute regime (that is, when the distance between the resonators is an order of magnitude greater than their size) we give explicit expressions for its value. Guided by this knowledge of how the infinite (bulk) chains behave, in Section 4 we design a finite chain of resonator dimers that has a topologically protected edge mode. This configuration takes inspiration from the bulk-boundary correspondence in the SSH model by introducing an interface, on either side of which the resonator dimers can be associated with different Zak phases thus creating a topologically protected edge mode.

In the quantum mechanical SSH model, the standard approach is to consider the tight-binding approximation. In this set-up, the Hamiltonian corresponding to the continuous differential problem is simplified by assuming that each particle only interacts with the surrounding crystal in a limited way that is easy to describe. This simplification gives a discrete approximation to the problem. Often, this is combined with a nearest-neighbour approximation, where long-range interactions are neglected, enabling explicit computations of the band structure. In this work, we prove that our system can be approximated by a discrete system, which captures all the interactions in full and has rigorous error estimates. In the dilute regime, we quantify the decay of the interactions and conclude that non-negligible interactions occur also for resonators separated by several unit cells. Since the edge modes are protected due to chiral symmetry, which is only present here under the nearest-neighbour approximation, we expect the midgap frequencies to be approximately stable with respect to errors which preserve this symmetry.

Finally, we conduct a fully-coupled numerical study of our finite chain of resonator dimers. This is based on an approach similar to that developed in [1]. We show that the crystal can exhibit topologically protected subwavelength edge modes in both the dilute and non-dilute regimes. Moreover, we study the stability of the midgap frequencies with respect to symmetry-preserving geometric errors. We show that, while the midgap frequency experiences variations (which is not the case under the nearest-neighbour approximation), these are much smaller than those seen by the band frequencies and the edge mode remains localized in the middle of the band gap even for very large geometric errors. We also make the comparison with a classical, unprotected, defect mode, similar to that studied in [3]. We show that the new subwavelength crystal exhibits a mode with a similar degree of localization but with greatly improved stability with respect to errors.

2 Preliminaries

In this section, we briefly review the layer potential operators and Floquet-Bloch theory that will be used in the subsequent analysis. More details on this material can be found in, for example, [4]. We also briefly review topological properties of the band structure.

2.1 Layer potential techniques

Let $D\in\mathbb{R}^{3}$ be a bounded, multiply-connected domain with $N$ simply-connected components $D_{i}$. Further, suppose that there exists some $0<s<1$ so that $\partial D_{i}$ is of class $C^{1,s}$ for each $i=1,\ldots,N$. Let $G^{0}$ and $G^{k}$ be the Laplace and (outgoing) Helmholtz Green’s functions, respectively, defined by

[TABLE]

We introduce the single layer potential $\mathcal{S}_{D}^{k}:L^{2}(\partial D)\rightarrow H_{\textrm{loc}}^{1}(\mathbb{R}^{3})$, defined by

[TABLE]

Here, the space $H_{\textrm{loc}}^{1}(\mathbb{R}^{3})$ consists of functions that are square integrable on every compact subset of $\mathbb{R}^{3}$ and have a weak first derivative that is also square integrable. It is well-known that $\mathcal{S}_{D}^{0}:L^{2}(\partial D)\rightarrow H^{1}(\partial D)$ is invertible, where $H^{1}(\partial D)$ is the set of functions that are square integrable on $\partial D$ and have a weak first derivative that is also square integrable.

We also define the Neumann-Poincaré operator $\mathcal{K}_{D}^{k,*}:L^{2}(\partial D)\rightarrow L^{2}(\partial D)$ by

[TABLE]

where $\partial/\partial\nu_{x}$ denotes the outward normal derivative at $x\in\partial D$.

The following relations, often known as jump relations, describe the behaviour of $\mathcal{S}_{D}^{k}$ on the boundary $\partial D$ (see, for example, [4]):

[TABLE]

and

[TABLE]

where $|_{\pm}$ denote the limits from outside and inside $D$.

2.2 Floquet-Bloch theory and quasiperiodic layer potentials

A function $f(x)\in L^{2}(\mathbb{R})$ is said to be $\alpha$-quasiperiodic, with quasiperiodicity $\alpha\in\mathbb{R}$, if $e^{-\mathrm{i}\mkern 1.0mu\alpha x}f(x)$ is periodic. If the period is $L>0$, the natural space for $\alpha$ is $Y^{*}:=\mathbb{R}/\tfrac{2\pi}{L}\mathbb{Z}\simeq(-\pi/L,\pi/L]$. $Y^{*}$ is known as the first Brillouin zone. Given a function $f\in L^{2}(\mathbb{R})$, the Floquet transform is defined as

[TABLE]

$\mathcal{F}[f]$ is always $\alpha$-quasiperiodic in $x$ and periodic in $\alpha$. Let $Y_{0}=[-L/2,L/2)$ be the unit cell. The Floquet transform is an invertible map $\mathcal{F}:L^{2}(\mathbb{R})\rightarrow L^{2}(Y_{0}\times Y^{*})$, with inverse (see, for instance, [4, 30])

[TABLE]

where $g(x,\alpha)$ is the quasiperiodic extension of $g$ for $x$ outside of the unit cell $Y_{0}$.

We will consider a three-dimensional problem which is periodic in one dimension. Define the unit cell $Y$ as $Y:=Y_{0}\times\mathbb{R}^{2}$. We define the quasiperiodic Green’s function $G^{\alpha,k}(x,y)$ as the Floquet transform of $G^{k}(x,y)$ in the first dimension of $x$, i.e.,

[TABLE]

Let $D$ be as in the previous layer potential definitions, but assume additionally $D\Subset Y$. We define the quasiperiodic single layer potential $\mathcal{S}_{D}^{\alpha,k}$ by

[TABLE]

It is known that $\mathcal{S}_{D}^{\alpha,0}:L^{2}(\partial D)\rightarrow H^{1}(\partial D)$ is invertible if $\alpha\neq 0$ [4]. It satisfies the jump relations

[TABLE]

and

[TABLE]

where $(\mathcal{K}_{D}^{-\alpha,k})^{*}$ is the quasiperiodic Neumann-Poincaré operator, given by

[TABLE]

2.3 Band structure and topological properties

In this section we briefly review the topological nature of the Bloch eigenbundle. For further details, and discussions of the topological quantities involved, we refer to [12, 14, 29]. Let $\mathcal{L}$ be a linear elliptic differential operator which is self-adjoint in $L^{2}(\mathbb{R}^{3})$ and whose coefficients are periodic in one dimension. Denote by $\mathcal{L}(\alpha)$ the operator with the same differential expression but acting on $\alpha$-quasiperiodic functions. It is well-known [31] that the spectrum $\sigma(\mathcal{L})$ of the original operator can be expressed in terms of the spectra $\sigma(\mathcal{L}(\alpha))$ as

[TABLE]

This describes a band structure of the spectrum of $\mathcal{L}$: for each $\alpha$ the spectrum $\sigma(\mathcal{L}(\alpha))$ is known to be discrete and will thus trace out bands $\sigma_{n}(\mathcal{L}(\alpha)),n=1,2,\dots$, as $\alpha$ varies. The spectrum of $\mathcal{L}$ is said to have a band gap if, for some $n$, $\max_{\alpha}\sigma_{n}(\mathcal{L}(\alpha))<\min_{\alpha}\sigma_{n+1}(\mathcal{L}(\alpha))$. A band is said to be non-degenerate if it does not intersect any other band.

On a non-degenerate band, indexed by $n=1,2,\dots$, there exists a family of associated Bloch eigenmodes $\{u_{n}^{\alpha}\}_{\alpha\in Y^{*}}$ which we define so that they are both normalized and depend continuously on $\alpha$. Observe that the base space $Y^{*}$ has the topology of a circle. A natural question to ask, when considering the topological properties of a crystal, is whether properties are preserved after parallel transport around $Y^{*}$. In particular, a powerful quantity to study is the Berry-Simon connection $A_{n}$, defined as

[TABLE]

For any $\alpha_{1},\alpha_{2}\in Y^{*}$, the parallel transport from $\alpha_{1}$ to $\alpha_{2}$ is $u_{n}^{\alpha_{1}}\mapsto e^{\mathrm{i}\mkern 1.0mu\theta}u_{n}^{\alpha_{2}}$, where $\theta$ is given by

[TABLE]

Thus, it is enlightening to introduce the so-called Zak phase, $\varphi_{n}^{z}$, defined as

[TABLE]

which corresponds to parallel transport around the whole of $Y^{*}$. When $\varphi_{n}^{z}$ takes a value that is not a multiple of $2\pi$, we see that the eigenmode has gained a non-zero phase after parallel transport around the circular domain $Y^{*}$. In this way, the Zak phase captures topological properties of the crystal. For crystals with inversion symmetry, the Zak phase is known to only attain the values [math] or $\pi$ [19, 47].

Remark 2.1**.**

In quantum mechanical contexts, $\langle\cdot,\cdot\rangle$ is typically chosen as the $L^{2}(Y)$-inner product. When working with Helmholtz scattering problems, however, this choice is not appropriate since the solutions are not elements of $L^{2}(Y)$. Instead, we will work with the $L^{2}(D)$-inner product. We will see that the behaviour on $D$ is enough to characterize non-trivial topological behaviour and capture the structure’s wave propagation properties.

3 Infinite, periodic chains of subwavelength resonator dimers

In this section, we study a periodic arrangement of subwavelength resonator dimers. This is an analogue of the SSH model. The goal is to derive a topological invariant which characterises the crystal’s wave propagation properties and indicates when it supports topologically protected edge modes. The analysis here holds for a very general class of high-contrast resonator chains, requiring only two assumptions of geometric symmetry.

3.1 Problem description

Assume we have a one-dimensional crystal in $\mathbb{R}^{3}$ with repeating unit cell $Y:=[-\frac{L}{2},\frac{L}{2}]\times\mathbb{R}^{2}$. Each unit cell contains two resonators (often referred to as a dimer) surrounded by some background medium. Suppose the resonators together occupy the domain $D:=D_{1}\cup D_{2}$. As well as sufficient smoothness for the above layer potential operators to be well defined, we need two assumptions of symmetry for the analysis that follows. The first is that each individual resonator is symmetric in the sense that there exists some $x_{1}\in\mathbb{R}$ such that

[TABLE]

where $R_{1}$ and $R_{2}$ are the reflections in the planes $p_{1}=\{-x_{1}\}\times\mathbb{R}^{2}$ and $p_{2}=\{x_{1}\}\times\mathbb{R}^{2}$, respectively. We also assume that the dimer is symmetric in the sense that

[TABLE]

Denote the full crystal by $\mathcal{C}$, that is,

[TABLE]

We denote the separation of the resonators within each unit cell, along the first coordinate axis, by $d:=2x_{1}$ and the separation across the boundary of the unit cell by $d^{\prime}:=L-d$.

Wave propagation inside the infinite periodic structure is modelled by the Helmholtz problem

[TABLE]

Here, $\omega$ is the frequency of the incident waves which is assumed to be small, such that we are in a subwavelength regime. We refer to [4] for the definition of the outgoing radiation condition. The material parameter $\delta$ represents the contrast between the resonators and the background. In order for subwavelength resonant modes to exist, we assume that $\delta$ satisfies the high-contrast condition

[TABLE]

As an example, in the case of acoustic wave propagation, $\delta=\rho_{r}/\rho_{0}$ is the density contrast between the resonator material and the background material.

Let $\sigma$ be the spectrum of the operator

[TABLE]

acting on functions which satisfy the outgoing radiation condition in $x_{2},x_{3}$. Here, $\chi_{\mathcal{C}}$ denotes the indicator function of the periodic crystal $\mathcal{C}$. Only frequencies $\omega$ such that $\omega^{2}\in\sigma$ can be solutions to (3.4). Any other frequencies are not able to propagate in the material. It is worth emphasizing that, due to radiation in $x_{2}$- and $x_{3}$-directions, the resonant frequencies are complex with negative imaginary parts. Nevertheless, as we will see in Theorem 3.4, the resonant frequencies are real at leading order so we consider only their real parts in this work.

By applying the Floquet transform, the Bloch eigenmode $u_{\alpha}(x):=\mathcal{F}[u](x,\alpha)$ is the solution to the Helmholtz problem

[TABLE]

We refer to [4] for the definition of the $\alpha$-quasiperiodic outgoing radiation condition. We denote by $\sigma_{\alpha}$ the spectrum of the operator

[TABLE]

acting on $\alpha$-quasiperiodic functions which satisfy the outgoing radiation condition in $x_{2},x_{3}$. As discussed in Section 2.3, we have

[TABLE]

which describes the band structure of the crystal.

3.2 Analysis of quasiperiodic problem

In this section we conduct a thorough analysis of the band structure and topological properties of (3.6). We use the methods from [5, 6] to formulate the quasiperiodic resonance problem as an integral equation. The solution $u_{\alpha}$ of (3.6) can be represented as

[TABLE]

for some density $\Psi^{\alpha}\in L^{2}(\partial D)$. Then, using the jump relations (2.4) and (2.5), it can be shown that (3.6) is equivalent to the boundary integral equation

[TABLE]

where

[TABLE]

3.2.1 Quasiperiodic capacitance matrix

Let $V_{j}^{\alpha}$ be the solution to

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. We then define the quasiperiodic capacitance matrix $C^{\alpha}=(C_{ij}^{\alpha})$ by

[TABLE]

We will see, shortly, that finding the eigenpairs of this matrix represents a leading order approximation to the differential problem (3.6). First, however, we show some useful properties of $C^{\alpha}$.

Lemma 3.1**.**

The matrix $C^{\alpha}$ is Hermitian with constant diagonal, i.e.,

[TABLE]

Proof.

From the definition (3.10), it clearly follows that $C^{\alpha}$ is Hermitian. To show that $C_{11}^{\alpha}=C_{22}^{\alpha}$, we define the mapping $T$ by

[TABLE]

Then, thanks to the assumed symmetry of the dimer (3.2), it holds that $TV_{1}^{\alpha}=\overline{V_{2}^{\alpha}}$ and $TV_{2}^{\alpha}=\overline{V_{1}^{\alpha}}$. It follows that

[TABLE]

Using the jump conditions, in the case $\alpha\neq 0$, it can be shown that the capacitance coefficients $C_{ij}^{\alpha}$ are also given by

[TABLE]

where $\psi_{j}^{\alpha}$ are defined by

[TABLE]

Since $C^{\alpha}$ is Hermitian, the following lemma follows directly.

Lemma 3.2**.**

The eigenvalues and corresponding eigenvectors of the quasiperiodic capacitance matrix are given by

[TABLE]

where, for $\alpha$ such that $C_{12}^{\alpha}\neq 0$, $\theta_{\alpha}\in[0,2\pi)$ is defined to be such that

[TABLE]

In the dilute regime, we are able to compute asymptotic expansions for the band structure and topological properties. In this regime, we assume that the resonators can be obtained by rescaling fixed domains $B_{1},B_{2}$ as follows:

[TABLE]

for some small parameter $\varepsilon>0$.

We introduce the capacitance $\textrm{Cap}_{B}$ of the fixed domains as follows. Let $B=B_{i}$ for $i=1$ or $i=2$ and define

[TABLE]

where $\phi_{B}:=(\mathcal{S}_{B}^{0})^{-1}[\chi_{\partial B}]$. Due to symmetry, the capacitance is the same for the two choices $i=1,2$. It is easy to see that, by a scaling argument,

[TABLE]

Lemma 3.3**.**

We assume that the resonators are in the dilute regime specified by (3.13). We also assume that $\alpha\neq 0$ is fixed. Then we have the following asymptotics of the capacitance matrix $C_{ij}^{\alpha}$ as $\varepsilon\rightarrow 0$:

[TABLE]

Taking the imaginary part of (3.16), the corresponding asymptotic formula holds uniformly in $\alpha\in Y^{*}$.

Proof.

Recall that the capacitance matrix $C_{ij}^{\alpha}$ can be written as

[TABLE]

for $\alpha\neq 0$. We shall compute the asymptotics of $(\mathcal{S}_{D}^{\alpha,0})^{-1}$ for small $\varepsilon$.

Let us decompose the Green’s function $G^{\alpha,k}$ as

[TABLE]

where

[TABLE]

Then let us define

[TABLE]

Note that ${S}_{D}^{\alpha,0}={S}_{D}^{0}+\widetilde{S}_{D}^{\alpha,0}$.

Let us write the quasiperiodic single layer potential $\mathcal{S}_{D}^{\alpha,0}$ in a matrix form as

[TABLE]

and then decompose it as

[TABLE]

In order to keep the order of the norms in $L^{2}(\partial D)$ and $H^{1}(\partial D)$ constant as $\varepsilon\rightarrow 0$, we let $\mathcal{L}$ and $\mathcal{H}$, respectively, denote the spaces $L^{2}(\partial D)$ and $H^{1}(\partial D)$ along with the inner products

[TABLE]

Then, for a fixed $\widetilde{\varphi}\in L^{2}(\partial B)$, if we define $\varphi\in L^{2}(\varepsilon\partial B)$ as $\varphi(x)=\widetilde{\varphi}(\varepsilon^{-1}x)$ we have $\|\varphi\|_{\mathcal{L}}=O(1)$ as $\varepsilon\rightarrow 0$.

Next, we estimate the operator norms of $S_{I}$ and $S_{II}$. We first handle the operator $S_{I}$. By the scaling property $\mathcal{S}_{\varepsilon B}^{0}[\varphi]=\varepsilon\mathcal{S}_{B}^{0}[\widetilde{\varphi}],$ it can be shown that

[TABLE]

which implies that

[TABLE]

Here, the notation $A\lesssim B$ means that there exists a constant $K$ independent of $\varepsilon$ such that $A\leq KB$ for all small enough $\varepsilon$. Further, $\mathcal{B}(X,Y)$ is used to denote the space of bounded linear operators between the normed spaces $(X,\|\cdot\|_{X})$ and $(Y,\|\cdot\|_{Y})$, and the $\|\cdot\|_{\mathcal{B}(X,Y)}$ norm is defined in the usual way as $\|T\|_{\mathcal{B}(X,Y)}:=\inf\{M:\|T(x)\|_{Y}\leq M\|x\|_{X},\forall x\in X\}$.

Let us now consider $S_{II}$. We introduce the notation

[TABLE]

We have, for small $\varepsilon$, that

[TABLE]

Here, $y_{m}^{*}$ means a point on the line segment joining $y$ and $z_{j}$. Note that the series in the remainder term converges. Moreover, the gradient of the remainder term is of the same order. Since $\int_{\partial D_{j}}\varphi\>\mathrm{d}\sigma=O(\varepsilon^{2}\|\varphi\|_{\mathcal{L}})$, we get

[TABLE]

Here, $\nabla_{\partial D}$ is used to denote the surface gradient on $\partial D$. From this we have

[TABLE]

Similarly, we can show that

[TABLE]

and \|{\mathcal{S}}_{D_{j}}^{0}\big{|}_{\partial D_{i}}\|_{\mathcal{B}(\mathcal{L},\mathcal{H})}\lesssim\varepsilon^{2}. These imply that

[TABLE]

We now compute the asymptotic behaviour of $(\mathcal{S}_{D}^{\alpha,0})^{-1}$. We use the definition $\phi_{j}:=S_{I}^{-1}[\chi_{\partial D_{j}}]$ and introduce the capacitance of each individual resonator $D_{j}$ as $\textrm{Cap}_{D_{j}}:=\nolinebreak-\int_{\partial D_{j}}\phi_{j}\>\mathrm{d}\sigma.$ Note that $\textrm{Cap}_{D_{j}}=\varepsilon\textrm{Cap}_{B}$ by (3.14). Since we know from (3.17) and (3.20) that $S_{I}^{-1}S_{II}=O(\varepsilon)$ in the operator norm, applying the Neumann series gives

[TABLE]

We also have

[TABLE]

Then, from (3.18) and (3.19), together with the fact that $\|\phi_{j}\|_{\mathcal{L}}=\|\mathcal{S}_{I}^{-1}[\chi_{\partial D_{j}}]\|_{\mathcal{L}}\lesssim\varepsilon^{-1},\ j=1,2$, we obtain the series representations

[TABLE]

We are ready to compute the capacitance matrix. From (3.21), we have

[TABLE]

The expression for $C_{12}^{\alpha}$ can be derived in the same way. ∎

3.2.2 Band structure and Bloch eigenmodes

Define normalized extensions of $V_{j}^{\alpha}$ as

[TABLE]

where $|D_{1}|$ is the volume of one of the resonators ($|D_{1}|=|D_{2}|$ thanks to the dimer’s symmetry (3.2)). Using similar arguments to those given in [5, 8, 10], the following two approximation results can be proved.

Theorem 3.4**.**

The characteristic values $\omega_{j}^{\alpha}=\omega_{j}^{\alpha}(\delta),\leavevmode\nobreak\ j=1,2$, of the operator $\mathcal{A}^{\alpha}(\omega,\delta)$, defined in (3.8), can be approximated as

[TABLE]

where $\lambda_{j}^{\alpha},\leavevmode\nobreak\ j=1,2$, are eigenvalues of the quasiperiodic capacitance matrix $C^{\alpha}$.

Theorem 3.5**.**

The Bloch eigenmodes $u_{j}^{\alpha},\leavevmode\nobreak\ j=1,2$, corresponding to the resonances $\omega_{j}^{\alpha}$, can be approximated as

[TABLE]

where $\left(\begin{smallmatrix}a_{j}\\ b_{j}\end{smallmatrix}\right),\leavevmode\nobreak\ j=1,2,$ are the eigenvectors of the quasiperiodic capacitance matrix $C^{\alpha}$, as given by 3.2.

Remark 3.6**.**

From Theorems 3.4 and 3.5, we can see that the capacitance matrix can be considered to be a discrete approximation of the differential problem (3.6), since its eigenpairs directly determine the resonant frequencies and the Bloch eigenmodes (at leading order). This is analogous to the tight-binding model commonly used in the quantum-mechanical SSH system.

Remark 3.7**.**

In the quantum-mechanical SSH model, the tight-binding model is typically handled with a nearest-neighbour approximation, where only the interactions between neighbouring particles are considered. In this regime, the model is described by the simple Hamiltonian matrix

[TABLE]

where $v$ and $w$ are the two inter-particle coupling constants. Compare this to our discrete approximation, given by the capacitance matrix $C^{\alpha}$. If we applied a nearest-neighbour approximation, the capacitance matrix $C^{\alpha}$ would have the same form as the Hamiltonian (3.22) (up to an additive diagonal matrix). This would be achieved by neglecting the series in (3.15) and truncating the series in (3.16) to $|m|\leq 1$ only. However, in classical wave propagation problems such as these, the slow decay of the capacitance matrix means this approximation may be inaccurate. This is because non-negligible interactions exist even between resonators separated by several unit cells. We shall see that this is also the case for finite crystals, in Section 4.2.1.

Remark 3.8**.**

Theorem 3.5 shows that the Bloch eigenmodes are asymptotically constant on each resonator. The value attained on each successive resonator differs by a phase factor determined by $\theta_{\alpha}$. This theorem was proved in [10] using layer potential techniques under the assumption $\alpha\neq 0$. By analytic continuation we may extend to any $\alpha\in Y^{*}$ [14].

We now introduce notation which, thanks to the assumed symmetry of the resonators, will allow us to prove topological properties of the chain. Divide $Y$ into two subsets $Y=Y_{1}\cup Y_{2}$, where $Y_{1}:=[-\frac{L}{2},0]\times\mathbb{R}^{2}$ and let $Y_{2}:=[0,\frac{L}{2}]\times\mathbb{R}^{2}$, as depicted in Figure 2. Define $q_{1}$ and $q_{2}$ to be the central planes of $Y_{1}$ and $Y_{2}$, that is, the planes $q_{1}:=\{-\frac{L}{4}\}\times\mathbb{R}^{2}$ and $q_{2}:=\{\frac{L}{4}\}\times\mathbb{R}^{2}$. Let $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ be reflections in the respective planes. Observe that, thanks to the assumed symmetry of each resonator (3.1), the “complementary” dimer $D^{\prime}=D_{1}^{\prime}\cup D_{2}^{\prime}$, given by swapping $d$ and $d^{\prime}$, satisfies $D_{i}^{\prime}=\mathcal{R}_{i}D_{i}$ for $i=1,2$. Define the operator $T_{\alpha}$ on the set of $\alpha$-quasiperiodic functions $f$ on $Y$ as

[TABLE]

where the factor $e^{-\mathrm{i}\mkern 1.0mu\alpha L}$ is chosen so that the image of a continuous ($\alpha$-quasiperiodic) function is continuous.

We will now proceed to use $T_{\alpha}$ to analyse the different topological properties of the two dimer configurations. Define the quantity ${C_{12}^{\alpha}}^{\prime}$ analogously to $C_{12}^{\alpha}$ but on the dimer $D^{\prime}$, that is, to be the top-right element of the corresponding quasiperiodic capacitance matrix, defined in (3.10).

Lemma 3.9**.**

We have

[TABLE]

Consequently, if $d=d^{\prime}=\frac{L}{2}$ then $C_{12}^{\pi/L}=0$.

Proof.

Define ${V_{1}^{\alpha}}^{\prime},{V_{2}^{\alpha}}^{\prime}$ by (3.9) but on $D^{\prime}$ instead of $D$. Observe that $T_{\alpha}{V_{1}^{\alpha}}^{\prime}=e^{-\mathrm{i}\mkern 1.0mu\alpha L}\overline{V_{1}^{\alpha}}$ and $T_{\alpha}{V_{2}^{\alpha}}^{\prime}=\overline{V_{1}^{\alpha}}$. Then, we find that

[TABLE]

At $\alpha=\pi/L$, we have ${C_{12}^{\pi/L}}^{\prime}=-C_{12}^{\pi/L}$. Moreover, if $d=d^{\prime}$, the symmetry of the structure means that ${C_{12}^{\pi/L}}^{\prime}=C_{12}^{\pi/L}$ so it must be the case that $C_{12}^{\pi/L}=0$. ∎

Lemma 3.10**.**

We assume that $D$ is in the dilute regime specified by (3.13). Then, for $\varepsilon$ small enough,

• (i)

$\mathrm{Im}\ C_{12}^{\alpha}>0$* for $0<\alpha<\pi/L$ and $\mathrm{Im}\ C_{12}^{\alpha}<0$ for $-\pi/L<\alpha<0$. In particular, $\mathrm{Im}\ {C_{12}^{\alpha}}$ is zero if and only if $\alpha\in\{0,\pi/L\}$.*

• (ii)

$C_{12}^{\alpha}$* is zero if and only if both $d=d^{\prime}$ and $\alpha=\pi/L$.*

• (iii)

$C_{12}^{\pi/L}<0$* when $d<d^{\prime}$ and $C_{12}^{\pi/L}>0$ when $d>d^{\prime}$. In both cases we have $C_{12}^{0}<0$.*

The proof of 3.10 is given in Appendix B. This lemma describes the crucial properties of the behaviour of the curve $\{C_{12}^{\alpha}:\alpha\in Y^{*}\}$ in the complex plane. The periodic nature of $Y^{*}$ means that this is a closed curve. Part (i) tells us that this curve crosses the real axis in precisely two points. Taken together with (iii), we know that this curve winds around the origin in the case $d>d^{\prime}$, but not in the case $d<d^{\prime}$.

Theorem 3.11**.**

If $d\neq d^{\prime}$ there exists a band gap, for $\alpha$ away from zero. That is, for any small $\alpha_{0}>0$, we have that

[TABLE]

for small enough $\varepsilon$ and $\delta$.

The proof of 3.11 is given in Appendix C. The argument is based on representing the first and second resonant frequencies as

[TABLE]

and making use of the fact that, in the dilute regime and for fixed $\alpha\neq 0$, the capacitance coefficients can be expanded using 3.3.

Remark 3.12**.**

Part (ii) of 3.10 is a particularly deep result which shows that the dilute crystal has a degeneracy precisely when $d=d^{\prime}$. The methods developed in [5] can be applied to show that the dispersion relation has a Dirac cone at $\alpha=\pi/L$ in this case. As $d$ increases across the point where $d=d^{\prime}$, the band gap closes (to form a Dirac cone) before reopening. In 3.14 we will show that the reopened band gap has a non-trivial topology, similar to what has been observed in other systems (for example, in [45]).

Combining the results of 3.10, 3.2 and 3.5, we obtain the following result concerning the band inversion that takes place between the two geometric regimes $d<d^{\prime}$ and $d>d^{\prime}$.

Proposition 3.13**.**

For $\varepsilon$ small enough, the band structure at $\alpha=\pi/L$ is inverted between the $d<d^{\prime}$ and $d>d^{\prime}$ regimes. In other words, the eigenfunctions associated with the first and second bands at $\alpha=\pi/L$ are given, respectively, by

[TABLE]

when $d<d^{\prime}$ and by

[TABLE]

*when $d>d^{\prime}$. *

The eigenmode $S_{1}^{\pi/L}(x)+S_{2}^{\pi/L}(x)$ is constant and attains the same value on both resonators, while the eigenmode $S_{1}^{\pi/L}(x)-S_{2}^{\pi/L}(x)$ has values of opposite sign on the two resonators. They therefore correspond, respectively, to monopole and dipole modes, and 3.13 shows that the monopole/dipole nature of the first two Bloch eigenmodes are swapped between the two regimes. We will now proceed to define a topological invariant which we will use to characterise the topology of a chain and prove how its value depends on the relative sizes of $d$ and $d^{\prime}$. This invariant is intimately connected with the band inversion phenomenon, and we will prove that it is non-trivial only if $d>d^{\prime}$.

Theorem 3.14**.**

We assume that $D$ is in the dilute regime specified by (3.13). Then the Zak phase $\varphi_{j}^{z},j=1,2$, defined by

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the $L^{2}(D)$-inner product, satisfies

[TABLE]

for $\varepsilon$ and $\delta$ small enough.

Proof.

We compute the Zak phase $\phi^{z}_{j},j=1,2,$ of the first and second band, respectively. Observe that

[TABLE]

and in $D$ we have

[TABLE]

for all $\alpha\in Y^{*}$. By Theorem 3.5 and Lemma 3.2, it follows that

[TABLE]

so the Zak phase is given by

[TABLE]

Since we know that $\varphi_{j}^{z}$ is an integer multiple of $\pi$, we have for small enough $\delta$ that

[TABLE]

This representation of $\varphi_{j}^{z}$, which is analogous to well-known results for Hamiltonian systems such as the SSH model [12], shows that the Zak phase is given by the change in the argument of $C_{12}^{\alpha}$ as $\alpha$ varies over the Brillouin zone $Y^{*}$. We can see from parts (i) and (iii) of 3.10 that the winding number of the origin depends on whether $d<d^{\prime}$ or $d>d^{\prime}$. In the two cases we have, respectively, $\left[\theta_{\alpha}\right]_{Y^{*}}=0$ and $\left[\theta_{\alpha}\right]_{Y^{*}}=-2\pi$. Therefore, if $\delta$ is small enough, we have that

[TABLE]

Remark 3.15**.**

The dilute assumption is not necessary to conclude that the Zak phase is non-trivial for certain configurations. Combining 3.9 and (3.23), both of which are valid without assumptions of diluteness, we find that ${\varphi_{j}^{z}}^{\prime}-\varphi_{j}^{z}=\pi$, where ${\varphi_{j}^{z}}^{\prime}$ is the Zak phase of the crystal with resonator separation $d^{\prime}$ instead of $d$ (as in Figure 2). The assumption of diluteness is invoked to prove part (ii) of 3.10, which shows that there are only two different topological regimes and that degeneracy occurs only at $d=d^{\prime}$. We conjecture that this is true even without the dilute assumption, in which case it is not hard to prove 3.13 and 3.14.

Theorem 3.14 shows that the Zak phase of the crystal is non-zero precisely when $d>d^{\prime}$. The bulk-boundary correspondence suggests that we can create topologically protected subwavelength edge modes by joining half-space subwavelength crystals, one with $\varphi_{j}^{z}=0$ and the other with $\varphi_{j}^{z}=\pi$. According to 3.15, this is also valid in the non-dilute case. In Section 4, we will study a finite chain that is designed with this principle in mind and demonstrate that it exhibits an edge mode that is stable with respect to symmetry-preserving imperfections.

3.3 Numerical computations

The band structure and the Bloch eigenmodes were computed using the multipole expansion method derived in Appendix A. This relies on the assumption that the resonators are spherical, which is a special case of the more general geometry considered above. As shown in Theorem 3.5, the Bloch eigenmodes are asymptotically constant on each domain $D_{i}$ and hence accurate and efficient computations can be achieved by approximating functions by only the first term in the multipole expansion.

All the numerical computations in this paper were performed for the example of acoustic waves being scattered by air bubbles in water. This is a classic example of subwavelength resonance, where the resonant frequency of a single bubble is known as the Minnaert resonance [2, 37]. Throughout this paper, we use $\delta=10^{-3}$, which is roughly the density contrast between air and water. We also use the material parameters $d=12,d^{\prime}=42,L=54$ to exemplify a dilute crystal, and parameters $d=3,d^{\prime}=6,L=9$ to exemplify a non-dilute crystal.

3.3.1 Band structure

Figures 4 and 4 show the band structure in a dilute and a non-dilute crystal, respectively. In the dilute case, the sound cone occupies a majority of the Brillouin zone, which is expected due to the lower interactions between the resonators. Both crystals show the opening of a band gap.

3.3.2 Bloch eigenmodes

Figure 5 shows the first two Bloch eigenmodes for the crystal at $\alpha=\pi/L$ in the cases $d<d^{\prime}$ and $d>d^{\prime}$. The band inversion is clearly seen: when $d>d^{\prime}$ the monopole/dipole modes correspond to the second/first mode, respectively. The band inversion property demonstrates the fact that the crystal has a non-zero Zak phase when $d>d^{\prime}$. As $\alpha$ varies, the phase shift $\theta_{\alpha}$ between the values of the eigenmodes winds around 0, resulting in band inversion at some point $\alpha\in Y^{*}$.

4 Finite chains of subwavelength resonators

In this section, we will study a finite chain of resonators which has been carefully designed to support topologically protected edge modes. Specifically, we assume that $D$ has the form

[TABLE]

where $D_{0}$ is a single repeating resonator. In other words, $D$ consists of an odd number $N$ of identical resonators ($N=4M+1$) with alternating distances $d$ and $d^{\prime}$ that are swapped at the middle resonator. An example of such a configuration is depicted in Figure 6. This structure is based on the intuition that if one joins together two chains with different topological properties, a protected edge mode will occur at the interface (this is the principle of bulk-boundary correspondence). In Figure 6 it is shown how on either side of the central resonator (which constitutes the “edge”) one can associate each successive pair of resonators with dimers belonging to infinite chains that have different Zak phases.

We model wave propagation in the crystal $D$ by the Helmholtz problem

[TABLE]

4.1 Integral equation formulation of the problem

The solution $u$ of (4.2) can be represented as

[TABLE]

for some density $\Psi\in L^{2}(\partial D)$. Then, analogous to the approach used in the quasiperiodic case in Section 3.2, the jump relations can be used to show that (4.2) is equivalent to the boundary integral equation

[TABLE]

where

[TABLE]

4.2 Capacitance matrix

Similar to the quasiperiodic case in Section 3.2.1, the resonant frequencies and eigenmodes of the finite chain can be expressed in terms of the capacitance matrix. Let $V_{j},j=1,...,N$ be the solution to

[TABLE]

We define the capacitance coefficients matrix $C=(C_{ij})$ by

[TABLE]

Once again, we can use the jump conditions to show that the capacitance coefficients $C_{ij}$ are also given by

[TABLE]

where the functions $\psi_{j}$ are defined by

[TABLE]

Observe that as $\delta\rightarrow 0$, we have $\lambda\rightarrow\tfrac{1}{2}^{+}$. Then, using Gohberg-Sigal theory for operator-valued functions [9, 25] we have the following lemma.

Lemma 4.1**.**

For any $\delta$ sufficiently small there are, up to multiplicity, $N$ characteristic values $\omega_{j}=\omega_{j}(\delta),j=1,...,N$, to the operator-valued analytic function $\mathcal{A}(\omega,\delta)$ such that $\omega_{j}(0)=0$ for all $j$ and $\omega_{j}$ depends on $\delta$ continuously.

The following theorem, proved in [7], shows that the eigenvalues of $C$ determine the resonance frequencies of the finite structure.

Theorem 4.2**.**

The characteristic values $\omega_{j}=\omega_{j}(\delta),\leavevmode\nobreak\ j=1,...,N$, of $\mathcal{A}(\omega,\delta)$ can be approximated as

[TABLE]

where $\lambda_{j},\leavevmode\nobreak\ j=1,...,N$, are eigenvalues of the capacitance matrix $C$ and $|D_{0}|$ is the volume of a single resonator.

4.2.1 Nearest-neighbour approximation

Drawing on parallels to how SSH chains are studied in quantum mechanics, an appealing approach to approximating the problem of wave scattering by a finite system of subwavelength resonators is to consider a nearest-neighbour approximation. That is, to disregard long-range interactions between resonators, instead only considering the interactions between neighbouring elements. Mathematically, this means approximating the capacitance matrix (4.5) by setting $C_{ij}=0$ if $|i-j|>1$, giving a tridiagonal matrix. Intuitively, one would expect that such an approach will only give a good estimate to the problem in the dilute regime.

We wish to prove estimates on the extent to which a tight-binding approach can approximate the problem in the dilute regime. We consider a dilute system by rescaling the canonical domain $D_{0}$ in (4.1) as $D_{0}=\varepsilon B$, where $B$ is some connected domain that has size of order one. In this dilute regime, we are able to obtain an explicit representation of the capacitance matrix $C$ for the finite system (4.1). As in Section 3.2.1, we denote the capacitance of the fixed domain $B$ by $\textrm{Cap}_{B}$.

Lemma 4.3**.**

Consider a dilute system of $N$ identical subwavelength resonators with size of order $\varepsilon$, given by

[TABLE]

where $0<\varepsilon\ll 1$ and $z_{j}$ represents the fixed position of each resonator. In the limit as $\varepsilon\rightarrow 0$ the asymptotic behaviour of the corresponding capacitance matrix is given by

[TABLE]

Proof.

The argument is very similar to that in 3.3. We first write the single layer potential $\mathcal{S}_{D}^{0}$ in a decomposed matrix form, as

[TABLE]

We can then proceed to use scaling arguments, as in 3.3, to find estimates for $\|S_{I}\|$, $\|S_{I}^{-1}\|$ and $\|S_{II}\|$, and then use these bounds to compute the asymptotic behaviour of $C$. ∎

Remark 4.4**.**

The explicit representations for $C_{ij}$ derived in 4.3, when used in the formula from 4.2, give approximations for the resonant frequencies in the dilute regime. Moreover, the associated eigenmodes can be approximated using the fact that the characteristic functions $\Psi_{j}$, defined for each $\omega_{j}$ in (4.3), satisfy

[TABLE]

where $(a_{j1},\dots,a_{jN})$ is the eigenvector of $C$ associated with the eigenvalue $\lambda_{j}$. This approach is particularly useful for performing efficient numerical computations.

One can see from Lemma 4.3 that, for $i\neq j$, $C_{ij}$ satisfies the slow decay property

[TABLE]

This indicates that, for a system of subwavelength resonators, the nearest-neighbour approximation may not give an accurate representation. This is a significant difference between the classical wave propagation problems studied here and the analogous applications of topological insulator theory in quantum mechanics, where nearest-neighbour approximations are commonplace.

4.2.2 Chiral symmetry and edge mode frequencies

A prominent topic in the discussion of the SSH model is the notion of chiral symmetry. This is a geometric property which a system is said to possess if there is an unitary matrix $\Sigma$ with $\Sigma^{2}=I$ such that the capacitance matrix $\widetilde{C}$ satisfies $\Sigma\widetilde{C}\Sigma=-\widetilde{C}$. The significance of this property is that a chirally symmetric matrix will have a symmetric spectrum. This is easily seen from the fact that if $(\lambda,v)$ is an eigenpair for $\widetilde{C}$ then so is $(-\lambda,\Sigma v)$. Finite chains that have an odd number of resonators (such as the example studied here, Figure 6) will have an odd number of resonant frequencies hence there must be a middle frequency. Thus, if one can design a chain which has a band gap (which we have a suggestion of how to do from Section 3) and is chirally symmetric, there must be a midgap frequency.

The reason we use the notation $\widetilde{C}$ for the capacitance matrix in this discussion is that in quantum mechanical settings it is customary to define the zero-energy state to be such that the diagonal entries of the Hamiltonian (which plays the analogous role of the capacitance matrix) vanish. Thus, one constructs a translated capacitance matrix $\widetilde{C}$ by subtracting the constant diagonal elements. For the crystal in Figure 6, we can use (4.6) to approximate $\widetilde{C}$ by a nearest-neighbour approximation: a bisymmetric, tridiagonal matrix with odd size and zero diagonal. Such a matrix is chirally symmetric, and therefore has a zero eigenvalue. This shows that the finite system has a midgap frequency, at leading order.

The key property of a topologically protected state is that it retains its properties when imperfections exist in the structure. In particular, a chirally symmetric structure will retain its chiral symmetry when errors are made in the position of the resonators. This is because such errors will not affect the diagonal entries of $\widetilde{C}$ and, away from the diagonals, $\widetilde{C}_{ij}$ and $\widetilde{C}_{ji}$ will experience the same effects. Since the nearest-neighbour approximation of the capacitance matrix is chirally symmetric, we expect the midgap frequencies to be approximately stable with respect to errors in resonator position.

In Figure 9(f) we show how the resonant frequencies given by a nearest-neighbour approximation to a dilute resonator chain vary when subjected to errors in the position of the resonators. We use the multipole expansion method outlined in Appendix A to calculate the capacitance matrix (4.5) then 4.2 to compute the resonant frequencies from its eigenvalues. The pertinent conclusion from this is that, under the nearest-neighbour approximation, the midgap frequency is perfectly stable (as predicted by the above discussion). This approximation should be compared to Figure 9(a), where the same simulations are performed on a fully-coupled chain. In light of the slow decay of the off-diagonal terms in the capacitance matrix (4.6), the differences between the behaviour of the approximated and fully-coupled models are unsurprising, even when simulations are performed in a very dilute regime.

4.3 Numerical illustrations

We now perform a series of numerical computations to illustrate the difference between the topologically protected subwavelength localized modes in the finite dimer chain (4.1) and conventional, unprotected, subwavelength localized modes. The unprotected mode we study is produced by taking an equally spaced chain of resonators and changing the radius of the central resonator, thus introducing a defect (often known as a point defect). This system, depicted in Figure 7, is the finite, one-dimensional equivalent of the system studied in [3], where the existence of a subwavelength localized mode was proved in the case of an infinite crystal.

As was the case for the infinite chain in Section 3.3, the following numerical results for the finite chains were calculated for the case of acoustic waves being scattered by (subwavelength) air bubbles in water. The details of discretizing the operator $\mathcal{A}(\omega,\delta)$ using the multipole expansion method are given in Appendix A.

4.3.1 Existence of localized modes

Figures 8(a), LABEL: and 8(b) show the localized modes for the dimer and point-defect chains respectively (whose geometries are depicted in Figures 6, LABEL: and 7). The configurations have been chosen to give roughly the same strength of the localization.

4.3.2 Stability with respect to errors

Finally, we study the stability of the edge mode frequency with respect to random, symmetry-preserving imperfections. In light of the discussion in Section 4.2.2, we add random errors to the positions of the resonators and repeatedly compute the resonant frequencies. In Figures 9(a), LABEL: and 9(b) we can see that, in both the dilute and non-dilute regimes, the structure supports a localized mode (depicted in Figure 8(a) for the dilute case) whose resonant frequency is in the middle of the band gap. In Table 9(c) it is demonstrated that in the two regimes the stability of each frequency with respect to the random errors is very similar in magnitude. The fact that the midgap frequency is consistently further from the edges of the band gap in the non-dilute case is merely a consequence of the gap being wider in this regime. In Figure 9(a) we present the same simulations for a very short dimer chain, with only nine resonators. We can see, once again, that there is a midgap frequency which is much more stable than the bulk frequencies.

Finally, we make a comparison with the conventional defect mode exhibited by the subwavelength point-defect chain (shown in Figure 7). It is clear from Figure 9(d) that, even for relatively small errors, the frequency associated with the point-defect mode exhibits poor stability and is easily lost amongst the bulk frequencies. The comparison between the robustness of the two designs is particularly eye-opening in light of the observation that the degree of wave localization is very similar. The new, dimerized design is equally capable of localizing waves at subwavelength scales but does so with spectacularly enhanced robustness.

5 Concluding remarks

In this work, we have, both analytically and numerically, studied a fully-coupled chain of subwavelength resonator dimers. We have shown that the infinite crystal exhibits a non-trivial Zak phase in certain resonator configurations. In the dilute regime, we have given explicit expressions for the Zak phase, proved the existence of a non-trivial band gap and shown that band inversion occurs between the two different phase regimes. Guided by these findings, we have designed a finite resonator chain that exhibits topologically protected edge modes at its centre. This was based on being able to associate the dimers on either side of this edge with different values of the Zak phase. We have shown numerically that the edge mode frequency is well-localized in the band gap and that, when errors are added to the positions of the resonators, the variance of this frequency is significantly lower than that of the bulk frequencies. Although much of the explicit analysis was performed on infinite chains, numerical experiments showed that our approach can be used to create topologically protected edge modes in structures that contain only very small numbers of resonators.

Appendix A Multipole expansion method in three dimensions

Here we derive the multipole expansion approximation of $\mathcal{S}_{D}^{k}$ and $\mathcal{S}_{D}^{\alpha,k}$ in three dimensions. The method is a generalization of the method in two dimensions given in Appendix C of [6]. The overarching principle is that when working on spherical domains, the action of the single layer potential on spherical basis functions has an explicit, analytic representation.

The goal is to discretize the equations (3.7) and (4.3). Observe that the operators $\mathcal{A}$ and $\mathcal{A}^{\alpha}$ can be written as

[TABLE]

and

[TABLE]

so it is enough to find a discretized representation of the single layer potentials $\mathcal{S}_{D}^{k}$ and $\mathcal{S}_{D}^{\alpha,k}$.

For a radially symmetric Helmholtz equation, it is well-known that the spherical waves $j_{l}(kr)Y_{l}^{m}(\theta,\phi)$ and $h_{l}^{(1)}(kr)Y_{l}^{m}(\theta,\phi)$ gives a basis of solutions in the polar coordinates $(r,\theta,\phi)$. Here $Y_{l}^{m}(\theta,\phi),l\in\mathbb{N},m=-l,...,l$, are the spherical harmonics and $j_{n}(kr),h_{n}^{(1)}(kr)$ are the spherical Bessel and Hankel functions of the first kind, respectively, defined by

[TABLE]

where $J_{n}$ and $H_{n}^{(1)}$ are the ordinary Bessel and Hankel functions of the first kind.

We begin by deriving the multipole expansion of the single layer potential $\mathcal{S}_{D}^{k}$. The spherical harmonics $Y_{l}^{m}$ form a basis of $L^{2}(\partial D)$, and we seek the expansion of $\mathcal{S}_{D}^{k}$ in this basis. Define $u:=\mathcal{S}_{D}^{k}[Y_{l}^{m}]$, which is the solution to

[TABLE]

The above equation can be easily solved by the separation of variables technique in polar coordinates. It gives

[TABLE]

where $c=-\mathrm{i}\mkern 1.0mukR^{2}$.

In order to handle problems posed on disjoint domains, we will need an addition theorem relating spherical waves centred around a translated origin to spherical waves around the original origin. Suppose we have a point with coordinates $x=(r,\theta,\phi)$ in the original system and $x^{\prime}=(r^{\prime},\theta^{\prime},\phi^{\prime})$ in the translated system, with the coordinate vectors related by $x=x^{\prime}+b$ for $b=(r_{b},\theta_{b},\phi_{b})$. Moreover, we assume $r^{\prime}<r_{b}$. Then the addition theorem reads [24]

[TABLE]

where the coefficients $A_{l^{\prime}m^{\prime}}^{lm}$ are given by

[TABLE]

Here, the coefficients $C(l,m,l^{\prime},m^{\prime},\lambda,\mu)$ are in turn given by

[TABLE]

where we denote by

[TABLE]

the Wigner $3j$ symbols. To simplify these expressions slightly, we assume that the original coordinate system is aligned such that $b$ points along the positive $z$-axis, i.e., $\theta_{b}=0$. In this case

[TABLE]

Substituting this into the expression for $A_{l^{\prime}m^{\prime}}^{lm}$ gives

[TABLE]

Now, we compute the quasiperiodic single layer potential $\mathcal{S}_{D}^{\alpha,k}[Y_{l}^{m}]$ in the case when $D$ consists of a single resonator. Since

[TABLE]

we have

[TABLE]

Here, $D+n$ means a translation of the disk $D$ by $(nL,0,0)$ and $(r_{n}^{\prime},\theta_{n}^{\prime},\phi_{n}^{\prime})$ are the spherical coordinates with respect to the centre of $D+n$.

Using the addition theorem (A.3) we have

[TABLE]

where $Q_{\lambda}^{\mu}$ is the one-dimensional lattice sum in three dimensions, defined by

[TABLE]

An efficient method for computing this lattice sum is given in [34].

We are now ready to compute $\mathcal{S}_{D}^{\alpha,k}$ in the case when $D\Subset Y$ consists of two resonators, centred at $(-x_{1},0,0)$ and $(x_{1},0,0)$, respectively. This is what we require in order to perform computations on the infinite chain in Section 3.3. By identifying $L^{2}(\partial D)=L^{2}(\partial D_{1})\times L^{2}(\partial D_{2})$ we have

[TABLE]

Here the operator \mathcal{S}_{D_{i}}^{\alpha,k}\big{|}_{\partial D_{j}}:L^{2}(\partial D_{i})\rightarrow L^{2}(\partial D_{j}),i,j=1,2 is the evaluation of $\mathcal{S}_{D_{i}}^{\alpha,k}$ on $\partial D_{j}$. To compute the multipole expansion of \mathcal{S}_{D_{1}}^{\alpha,k}\big{|}_{\partial D_{2}}, we again use the addition theorem. We have

[TABLE]

In order to simulate the finite resonator chain in Section 4.3, we must now perform similar computations for the operator $\mathcal{S}_{D}^{k}$ in the case when $D$ consists of $N$ resonators. We assume the resonators to be arranged collinearly along the $x_{1}$-axis. By identifying $L^{2}(\partial D)=L^{2}(\partial D_{1})\times\ldots\times L^{2}(\partial D_{N})$ we have

[TABLE]

where, as in the quasiperiodic case, \mathcal{S}_{D_{i}}^{k}\big{|}_{\partial D_{j}}:L^{2}(\partial D_{i})\rightarrow L^{2}(\partial D_{j}) is the evaluation of $\mathcal{S}_{D_{i}}^{k}$ on $\partial D_{j}$. This relies on the addition theorem once again. The diagonal terms are easily evaluated using (A.2) directly. Away from the diagonals, the addition theorem (A.3) gives that

[TABLE]

Appendix B Proof of 3.10

Part (i): $\text{Im }C_{12}^{\alpha}>0$ for $0<\alpha<\pi/L$ and $\text{Im }C_{12}^{\alpha}<0$ for $-\pi/L<\alpha<0$. In particular, $\text{Im }{C_{12}^{\alpha}}$ is zero if and only if $\alpha\in\{0,\pi/L\}$.

Recall, from 3.3, the following expansion of ${C_{12}^{\pi/L}}$ for fixed $\alpha\neq 0$ in the dilute regime:

[TABLE]

where, in order to simplify notation, we have defined the function $f$ as

[TABLE]

The imaginary part

[TABLE]

converges for all $\alpha\in Y^{*}$, and (B.1) is valid for imaginary parts also at $\alpha=0$.

We will express $f$ in terms of Lerch’s Transcendent function $\Phi(z,s,a)$, after having first reviewed some basic properties. For details we refer to [18]. $\Phi(z,s,a)$ is defined by the power series

[TABLE]

for $z\in\mathbb{C}$ where this series converges, and extended by analytic continuation. If $\text{Re}(s)>0,\text{Re}(a)>0$ and $z\in\mathbb{C}\setminus[1,\infty)$ we have the integral representation

[TABLE]

where $\Gamma$ is the Gamma function.

Now, from the definition of $f$ in (B.2), we have

[TABLE]

Using the integral representation (B.4) we get, after simplifications,

[TABLE]

The imaginary part satisfies

[TABLE]

At the points $\alpha=0$ and $\alpha=\pi/L$, the functions $V_{1}^{\alpha}$ and $V_{2}^{\alpha}$ are real-valued and hence $\text{Im }C_{12}^{\alpha}=0$. We will show that, for $\varepsilon$ small enough, this imaginary part is zero precisely when $\alpha\in\{0,\pi/L\}$. The integrand in (B.6) is positive, and hence $\text{Im }f(\alpha,d)=0$ if and only if $\alpha\in\{0,\pi/L\}$ for $\alpha\in Y^{*}=(-\pi/L,\pi/L]$. This shows that the leading order term of $\text{Im }C_{12}^{\alpha}$ is zero precisely when $\alpha\in\{0,\pi/L\}$. Moreover, we can easily verify from (B.6) that

[TABLE]

This shows that for small enough $\varepsilon$, the function $\text{Im }C_{12}^{\alpha}$ will be monotonic around $\alpha=0$ and $\alpha=\pi/L$. Since we know that $\alpha\in\{0,\pi/L\}$ are exact zeros of $\text{Im }C_{12}^{\alpha},$ these zeros will be isolated for small enough $\varepsilon$. It follows that, for small enough $\varepsilon$, $\text{Im }C_{12}^{\alpha}$ is zero precisely when $\alpha\in\{0,\pi/L\}$. Then, from (B.6) it follows that $\text{Im }C_{12}^{\alpha}>0$ for $0<\alpha<\pi/L$ and $\text{Im }C_{12}^{\alpha}<0$ for $-\pi/L<\alpha<0$.

Part (ii): $C_{12}^{\alpha}$ is zero if and only if both $d=d^{\prime}$ and $\alpha=\pi/L$.

By part (i), any zero must satisfy $\alpha=0$ or $\alpha=\pi/L$. We begin by excluding the case $\alpha=0$. As $\alpha\rightarrow 0$, it is known that the quasiperiodic capacitance of a single particle vanishes [4, 6]. In other words, we have, for the total capacitance of $D$,

[TABLE]

where the last equality follows since $C_{12}^{0}$ is real. Since $C_{11}^{\alpha}>0$ we have $C_{12}^{0}<0$.

We now turn to the case $\alpha=\pi/L$. We already know from 3.9 that ${C_{12}^{\pi/L}}=0$ if $d=d^{\prime}$. To show that this is the only zero of ${C_{12}^{\pi/L}}$, we will show that ${C_{12}^{\pi/L}}$ is strictly monotonic as a function of $d<L$. From the definition of $f$ we compute the derivative

[TABLE]

The sum is an alternating series, with decreasing terms and negative first term. Hence it converges to a negative value, and by (B.1) we have

[TABLE]

for $\varepsilon$ small enough. This shows that ${C_{12}^{\pi/L}}$ has a unique zero when $d=d^{\prime}$, which completes the proof of part (ii).

Part (iii): $C_{12}^{\pi/L}<0$ when $d<d^{\prime}$ and $C_{12}^{\pi/L}>0$ when $d>d^{\prime}$. In both cases we have $C_{12}^{0}<0$.

We already know, from the proof of part (ii), that ${C_{12}^{0}}<0$. The other conclusions, namely that ${C_{12}^{\pi/L}}<0$ if $d<d^{\prime}$ and ${C_{12}^{\pi/L}}>0$ if $d>d^{\prime}$, follow directly from (B.7). ∎

Appendix C Proof of 3.11

We begin by proving the following lemma.

Lemma C.1**.**

For every $a,b$ such that $-1\leq a<1$ and $0<b<1$, the following holds:

[TABLE]

Proof.

We will split into the cases $a>0$ and $a\leq 0$, and begin with $a>0$. The right-hand side can be written as follows:

[TABLE]

Indeed, the integrand has a primitive function

[TABLE]

which shows (C.2). Then we have, for $a>0$

[TABLE]

where the last step follows because $\sinh(bt)-\frac{b\sinh(t)}{\frac{1-a}{2}\cosh(t)+\frac{1+a}{2}}>0$ for all $t>0$ in the case $a>0$. This proves the lemma in the case $a>0$. We now turn to the case $a\leq 0$. It is easy to see that for every $b$,

[TABLE]

Moreover, we have

[TABLE]

where we have used a known value for the integral (for example found in [52]). Together with (C.3), this proves the lemma in the case $a\leq 0$. ∎

Proof of Theorem 3.11. We will show that there is a frequency $\omega_{0}$ such that

[TABLE]

For sufficiently small $\delta$, by 3.4 and 3.2, it is enough to show that there is a constant $C_{0}$ such that

[TABLE]

Define $C_{0}$ as

[TABLE]

that is, $C_{0}$ is defined as the leading order of the eigenvalues of $C^{\alpha}$ at the degenerate point $d=d^{\prime},\alpha=\pi/L$. The sum appearing in the expansion of $C_{11}^{\alpha}$ can be explicitly computed as

[TABLE]

Then we have

[TABLE]

for $\varepsilon$ small enough. This shows that, for sufficiently small $\varepsilon$,

[TABLE]

We now turn to the second inequality of (C.4). By (B.1) and (B.5) we have

[TABLE]

Recall that $0<d<L$, so we can apply (C.1) with $b=d/L$ or with $b=(1-d)/L$ and with $a=\cos(\alpha L)$. Expanding the absolute value and applying (C.1), we find after simplifications that

[TABLE]

Together with (C.5), this shows that

[TABLE]

for $\varepsilon$ small enough. We have thus proved (C.4), from which the theorem follows. ∎