{'fraction_non_alphanumeric': 0.04285882607073425, 'fraction_numerical': 0.015817999725075114, 'mean_word_length': 4.110492247478549, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 46, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Functional analytic methods for discrete approximations of subwavelength resonator systems

We survey functional analytic methods for studying subwavelength resonator systems. In particular, rigorous discrete approximations of Helmholtz scattering problems are derived in an asymptotic subwavelength regime. This is achieved by re-framing the Helmholtz equation as a non-linear eigenvalue problem in terms of integral operators. In the subwavelength limit, resonant states are described by the eigenstates of the generalized capacitance matrix, which appears by perturbing the elements of the kernel of the limiting operator. Using this formulation, we are able to describe subwavelength resonance and related phenomena. In particular, we demonstrate large-scale effective parameters with exotic values. We also show that these systems can exhibit localized and guided waves on very small length scales. Using the concept of topologically protected edge modes, such localization can be made robust against structural imperfections.


Wave manipulation at subwavelength scales
A widespread ambition in wave physics is to be able to manipulate waves at scales that are much smaller than their wavelengths. On the other hand, an intuitive physical paradigm is that the propagation of a wave is not significantly affected by small objects or inhomogeneities. In particular, if an object is much smaller than the incident wavelength, then it will typically have a negligible scattering effect. This simplified phenomenon is closely related to Abbe's famous diffraction limit, which describes how the resolution of imaging systems depends on the operating wavelength. In order to overcome this limit, and be able to manipulate waves at subwavelength scales, there is widespread interest in settings where small objects exhibit subwavelength resonance and strongly scatter waves with comparatively large wavelengths.
The first high-profile example of subwavelength resonance came in the setting of acoustics when Marcel Minnaert observed the resonance of small air bubbles in water [57]. The very large contrast between the material parameters of air and water is understood to be the crucial mechanism here [13,55]. This phenomenon has since been observed in a variety of other settings, such as Helmholtz resonators [27], plasmonic nanoparticles [25,26] and high-contrast dielectric particles [3,24]. In general, bounded material inclusions whose parameters differ greatly from the background medium and which experience subwavelength resonance will be referred to as subwavelength resonators in this work.
The value of subwavelength resonators is that they can be used as the building blocks for large, complex structures which can exhibit a variety of exotic and useful properties. These micro-structured materials are examples of metamaterials: materials with a repeating micro-structure that exhibit properties surpassing those of the individual building blocks [48]. The widespread interest in metamaterials began with the realisation that they could be designed to have effectively negative material parameters [60] and, as a result, could be used to design perfect lenses [63,58] as well as cloaking and shielding devices [56,1]. Moreover, due to the subwavelength nature of the resonance, these structures enable wave control on very small length scales. Most notably, waves can be confined or guided using very small devices [15,51,67].
More recently, the study of micro-structured resonant media has focussed on designing structures whose properties are robust with respect to imperfections in their construction. This is important for realising the applications of this theory, since small errors will be introduced during the manufacturing process and many of the properties of these micro-structured media are very sensitive. Developments in this area have been based on studying the topological properties of periodic structures to create so-called topologically protected modes [10]. These concepts have previously been widely studied in a variety of settings, most notably in quantum mechanics for the Schrödinger operator [40,41,38].

Analysis of scattering problems
There is a large body of work dedicated to studying the scattering of waves by a collection of objects [53]. A popular simplifying assumption is to consider scattering by circular or spherical inclusions. In the case of a single inclusion, characterisations of the scattered field can be obtained through the use of expansions in terms of Bessel functions or spherical harmonics [33,44]. Likewise, in the case of two spheres a bispherical coordinate system can be used to give explicit representations of solutions [12].
In order to study scatterers with a more general class of shapes, integral equation methods are commonly used [34]. Boundary integral formulations can be used to reduce the dimension of the scattering problem, by rephrasing it as a problem posed on the boundaries of the scatterers [22,13]. Similarly, approaches that use Lippmann-Schwinger representations to express solutions in terms of volume integrals have been used for both scalar models [55,3] and for the Maxwell equations [35,36]. The fundamental idea here is that by representing solutions using appropriate integral operators a scattering problem can be equivalently phrased as a non-linear eigenvalue problem. With this formulation, scattering resonances can be characterised as the poles of meromorphic operator-valued functions [39]. In some settings, this can be paired with a scattering matrix [11] or transfer matrix [52] formulation to give a concise description of the response of the system.
Given the multi-scale nature of subwavelength metamaterials, asymptotic techniques are often used to understand their properties. In particular, a common approach for studying subwavelength problems (a) A system of finitely many resonators. · · · · · · .
(b) A periodic system of resonators. Figure 1: The functional analytic method developed here is useful for studying scattering by a system of material inclusions, which act as subwavelength resonators in an appropriate high-contrast regime. We are able to derive concise asymptotic results in terms of the capacitance matrix for the case of either finitely many resonators or a periodically repeating array of finitely many resonators.
is to assume that the resonator is asymptotically small while the other material parameters are fixed [55,3]. This is convenient because it can be implemented easily via a change of variables to give a concise description of a structure that is significantly (in an asymptotic sense) smaller than the operating wavelength. Related to this, homogenization techniques are often used to describe effective properties of micro-structured media. However, standard homogenization techniques do not apply here since these phenomena are based on local resonance of the small repeating units [28,23]. One downside to modelling subwavelength resonators as being asymptotically small is that it can tend to simplify the otherwise exotic behaviour as it reduces the underlying mechanism to just a rescaling of the model. Conversely, in this work we instead fix the resonators' size and position and consider an asymptotic limit in the material contrast parameter (which, in the case of acoustic waves, describes the ratio of the density inside and outside the resonators). This has the fundamental difference that the limiting problem is not trivial and has a spectrum of eigenvalues that can be understood. The asymptotic perturbation theory of Gohberg and Sigal [42,16] can then be used to prove the existence of subwavelength resonant frequencies, which are defined as resonant frequencies which satisfy a given asymptotic condition. This approach reveals the fundamental differences between a system's subwavelength resonant modes and the higher-frequency resonances.

Functional analytic approach
We represent the subwavelength resonators as material inclusions D i for i ∈ I, where I ⊂ N is some index set. The number of connected components that make up D should either be finite or the geometry should be periodic, such that it is given by an array of finitely many resonators that repeats indefinitely, as illustrated in Figure 1. We will study scalar wave settings where we throughout use ω to denote the frequency of the waves. We are interested in understanding solutions to Helmholtz resonance problems of the form satisfies an outgoing radiation condition, (1.1) where D = ∪ i∈I D i and the dimension d ∈ {2, 3}. Such Helmholtz equations, which can be used to model acoustic and polarized electromagnetic waves, represent the simplest model for wave propagation that still exhibits the rich phenomena associated to subwavelength physics. We wish to characterise solutions to (1.1) in terms of the system's subwavelength resonant modes. The parameters δ i in (1.1) are of crucial importance, and can be interpreted as the material contrast. In order to achieve subwavelength resonance we will assume that these parameters are small, corresponding to a large contrast between the materials. So that we can perform concise asymptotics in terms of the material contrast, we will introduce the real-valued parameter δ := |δ 1 | and assume that δ i = O(δ) as δ → 0 for all i ∈ I. We will then make a definition of a resonant mode being subwavelength as an asymptotic property in terms of δ. Definition 1.1 (Subwavelength resonant frequency). Given δ > 0, a subwavelength resonant frequency ω = ω(δ) ∈ C is defined to be such that (i) there exists a non-trivial solution to (1.1), known as an associated resonant mode; (ii) ω depends continuously on δ and satisfies ω → 0 as δ → 0.
The starting point for using functional analytic methods to understand resonance problems is to reframe the problem (1.1) as an operator equation. For instance, we will show that finding a solution to (1.1) is equivalent to finding a (non-trivial) function Φ ∈ L 2 (∂D) such that an integral equation of the form is satisfied. In the example studied in this work, A(ω, δ) will be an operator L 2 (∂D) × L 2 (∂D) → H 1 (∂D) × L 2 (∂D) and Φ an element of L 2 (∂D) × L 2 (∂D). Here, H 1 is the usual Sobolev space of square-integrable functions whose weak derivative is square integrable. If X and Y are two Banach spaces, then we write L(X, Y ) to denote the space of bounded linear operators from X into Y . In this work, we are interested in the case that X and Y are themselves spaces of functions and we have the following definition to describe the zeros of an operator-valued function that maps into L(X, Y ).

Definition 1.2 (Characteristic value).
A point z 0 ∈ C is said to be a characteristic value of T : C → L(X, Y ), which is an operator-valued function of a complex variable, if there exists some φ ∈ X such that φ(z 0 ) = 0 and T (z 0 )φ(z 0 ) = 0.
Comparing Definitions 1.1 and 1.2, we see that finding a subwavelength resonant frequency of the system is equivalent to finding, for a given δ, a characteristic value ω of A(ω, δ) which is such that ω(δ) → 0 as δ → 0. Our approach to finding such solutions is to consider perturbations of the elements of the kernel of A(0, 0). We will see that this space has dimension equal to the number of distinct resonators in the structure. Once we understand ker(A(0, 0)), we can characterise characteristic values of A(ω, δ) for small ω and δ as perturbations of this space. This analysis is based on the asymptotic perturbation theory of Gohberg and Sigal [42,16] and allows us both to prove the existence of subwavelength resonant modes (satisfying Definition 1.1) and to derive asymptotic formulas for their value.
This functional analytic approach has been used to describe subwavelength resonance in a variety of different physical settings. For instance, it was used to characterise a system of subwavelength Helmholtz resonators in [27], plasmonic particles in [25] and high-contrast dielectric resonators in [3,24]. In this work, we will explore its use to study scattering by a high contrast material inclusion, such as an bubble in water for the case of acoustic waves. This approach was first developed in this setting by [13] but, as we shall see, has since been developed to cover a variety of different settings and applications.

Capacitance coefficients
In the high-contrast Helmholtz setting that we will consider here, the functional analytic method described above will yield an approximation in terms of capacitance coefficients. Capacitance coefficients have a long history in electrostatics, where they govern the relationship between the distributions of potential and charge in a system of conductors. In particular, Maxwell introduced the matrix of capacitance coefficients C ∈ R N ×N to be such that if V ∈ R N is the vector of potentials on a system of N conductors then Q = CV is the vector of charges on the conductors [54,37].
Capacitance coefficients appear in the setting of subwavelength Helmholtz problems when we describe the principal part of the meromorphic operator (A(ω, 0)) −1 , which is a finite-rank operator governing the perturbation of the kernel of A(0, 0). Using a pole-pencil decomposition, we are able to project the problem onto ker(A(0, 0)), which gives a finite-dimensional characterization in terms of the generalized capacitance matrix. These ideas are elaborated in Appendix A, and shows that the capacitance coefficients appear naturally from the functional analytic approach described in Section 1.3.
In this article, we will survey how the generalized capacitance matrix offers a rigorous and intuitive discrete approximation to subwavelength Helmholtz scattering and resonance problems. This gives v δ 7 v 7 leading-order asymptotic expressions for both resonant modes and scattered solutions in terms of the eigenvalues and eigenvectors of the generalized capacitance matrix, which are accompanied by precise error bounds. We will see that a wide variety of different applications and phenomena can be studied using the capacitance approximation, demonstrating the power of reducing a differential problem to a matrix approximation in this way.

Finite systems
In this section, we apply the functional analytic method set out in Section 1.3 to a system of finitely many resonators. In particular, we suppose that d ∈ {2, 3}, N ∈ N and D 1 , D 2 , . . . , D N ⊂ R d are disjoint, connected sets with boundaries in C 1,s for some 0 < s < 1. An example of the setting is sketched in Figure 2. We use v i to denote the wave speed in resonator D i , then k i = ω/v i is the wave number in D i , where ω is the frequency of the incoming wave. Similarly, the wave speed and wave number in the background medium are denoted by v and k. Then, we consider a Helmholtz resonance problem for the domain D = D 1 ∪ · · · ∪ D N : u(x) satisfies the Sommerfeld radiation condition, where the Sommerfeld radiation condition says that and guarantees that energy is radiated outwards by the scattered solution. We assume that all contrast parameters are small while the wave speeds have order 1. In other words, the parameter δ > 0 is such that In order to concisely represent the different δ i , we introduce the function δ(x) ∈ L 2 (∂D) as

Main results of the capacitance formulation
The main tool that allows us to reveal the resonant properties of the system D = D 1 ∪ · · · ∪ D N is the Helmholtz single layer potential. This is an operator S ω D : where H 1 loc (R d ) is the set of functions in H 1 (A) for all compact subsets A ⊂ R d , while G ω is the Helmholtz Green's function, given by Here, H (1) 0 is the Hankel function of the first kind and order zero. In the case ω = 0, G 0 is the Laplace Green's function given by The single layer potential is useful because it allows us to seek solutions to (2.1) of the form where ψ, φ ∈ L 2 (∂D) are density functions that need to be found. The value of this representation is that a solution of the form (2.8) necessarily satisfies the Helmholtz equations and the radiation condition in problem (2.1). Also, the different wave numbers k i inside D i have been taken into account by using different single layer potentials S ki D . We can collectively represent these single layer potentials through a single operator S ω D , defined piecewise for i = 1, ..., N as We emphasize that S 0 D = S 0 D . It now remains only to find ψ and φ such that the transmission conditions on ∂D are fulfilled. This can be achieved through the introduction of an additional integral operator, the Neumann-Poincaré operator associated to D. This is an operator K ω, * D on L 2 (∂D) which is defined as We define K ω, * D in the same spirit as S ω D , namely We are now able to describe how S ω D and its normal derivative behave on ∂D. In particular, it holds that for any ϕ ∈ L 2 (∂D) [22] S ω and where the subscripts + and − denote taking the limit from outside and inside the boundary ∂D, respectively. With the so-called jump conditions (2.12) and (2.13) in hand, we can derive the following lemma, which characterises the resonance problem (2.1) as a boundary integral equation.
Lemma 2.1. The Helmholtz problem (2.1) is equivalent to finding ψ, φ ∈ L 2 (∂D) such that where the operator A(ω, δ) : The approach outlined in Section 1.3 can now be applied to Lemma 2.1 to prove the existence of subwavelength resonances, as defined in Definition 1.1, and derive their asymptotic behaviour as δ → 0. The idea here is to study the kernel of A(0, 0), where and then understand how ker A(0, 0) is perturbed when δ and ω are non-zero. The following lemma describes the two operators that appear in A(0, 0), as given in (2.15).
From Lemma 2.2 we can see that A(0, 0) has an N -dimensional kernel. Therefore, ω = 0 is a characteristic value of A(ω, 0). Due to symmetry, the multiplicity of ω = 0 is, in fact, 2N . When the material parameters are real, it is easy to see that A(ω, δ) = A(−ω, δ), from which we can see that the resonant frequencies will be symmetric with respect to the imaginary axis, in the sense described in Lemma 2.3 (cf. the analysis of [39]). Lemma 2.3. The set of resonant frequencies is symmetric in the imaginary axis. In particular, if δ i , v i ∈ R for all i = 1, ..., N , and if ω is such that (2.14) is satisfied for some non-zero ψ, φ ∈ L 2 (∂D), then it will also hold that With Lemma 2.3 in mind, we will subsequently state results only for the resonant frequencies with non-negative real parts. We can now show the following two theorems, using the approach described in Section 1.3.
Definition 2.5 (Capacitance matrix). For a system of N ∈ N resonators D 1 , . . . , D N in R 3 we can define the capacitance matrix C = (C ij ) ∈ R N ×N to be the square matrix given by Due to the different material parameters inside each resonator we introduce the generalized capacitance matrix, which is the main quantity we use in order to describe the subwavelength resonators. Definition 2.6 (Generalized capacitance matrix). For a system of N ∈ N resonators D 1 , . . . , D N in R 3 we can define the generalized capacitance matrix, denoted by C = (C ij ) ∈ C N ×N , to be the square matrix given by Theorem 2.7. Let d = 3. Consider a system of N subwavelength resonators in R 3 . As δ → 0, the N subwavelength resonant frequencies satisfy the asymptotic formula where {λ n : n = 1, . . . , N } are the eigenvalues of the generalized capacitance matrix C ∈ C N ×N , which satisfy λ n = O(δ) as δ → 0.
Remark 2.8. The assumption that the dimension d = 3 in Definitions 2.5 and 2.6 is important as the Laplace single layer potential S 0 D is known to be invertible in this case. As we will see in Section 2.5, this is not generally the case when d = 2 meaning that the corresponding version of Theorem 2.7 is slightly less elegant. Corollary 2.9. Let d = 3. Let v n be the normalized eigenvector of C associated to the eigenvalue λ n . Then the normalized resonant mode u n associated to the resonant frequency ω n is given, as δ → 0, by The capacitance matrix is defined solely in terms of the kernel of the integral operators A(0, 0) and A * (0, 0). In fact, the adjoint K 0 D of the Neumann-Poincaré operator satisfies ker(− 1 where, as in Lemma 2.2, {ψ 1 , ..., ψ N } is a basis for ker(− 1 2 I + K 0, * D ). In Appendix A, we use these ideas to define the generalized capacitance matrix purely in terms of the integral operator A, thus providing a general method to study subwavelength resonance systems.

Properties of the capacitance matrix
Through Theorem 2.7 and Corollary 2.9, we have reduced the resonance problem (2.1) to a matrix eigenproblem for the generalized capacitance matrix C. We now wish to understand the properties of C.
where δ ij is the Kronecker delta. Then, the capacitance coefficients, defined in Definition 2.5, are given by From Lemma 2.11 emerges a slightly different explanation for why the capacitance matrix approximation works. In the limiting case, when δ i = 0 for all i, the Helmholtz problem (2.1) is reduced to a Neumann boundary value problem in the interior of D and a Dirichlet boundary value problem in the exterior of D. When ω = 0, the interior Neumann problem is solved by constant functions, meaning that the solution of the exterior Dirichlet boundary value problem is a linear combination of the functions V 1 , . . . , V N . The projection onto this finite dimensional space yields a leading-order approximation of the solution as ω, δ → 0, in the form of an eigenvalue problem for the generalized capacitance matrix. Lemma 2.11 is also useful as it allows us to immediately see, among other things, the symmetry of the capacitance matrix.  Figure 3: A system of two spherical resonators can be described using bispherical coordinates. Such a coordinate system is convenient since the boundaries of the spheres lie on level sets and the capacitance coefficients can be calculated explicitly.
Lemma 2.12. The capacitance matrix C is symmetric and positive definite.
The symmetry and positive definiteness of the capacitance matrix C is useful for understanding the properties of the generalised capacitance matrix C which is the product of C with a diagonal matrix containing the weights δ i v 2 i /|D i |. In the case that δ i v 2 i > 0, this diagonal matrix is positive definite so, for example, we have the following lemma, which can be proved using the fact that C is always Hermitian.
. . , N , are real-valued positive numbers, then the generalized capacitance matrix C ∈ C N ×N has N linearly independent eigenvectors. Remark 2.14. We will see, in Section 4.3, that non-zero imaginary parts of δ i v 2 i can be used to model damping and amplification in the system. In this case, we can create exceptional points where eigenvalues and eigenvectors coincide and C is not diagonalizable.
Explicit formulas for capacitance coefficients are generally beyond reach. However, by making some additional assumptions, we can make the capacitance matrix easier to understand. For example, spherical resonators can be easily described using appropriate radial coordinate systems. In this way, we can see that if D is a single sphere of radius r then its capacitance is given by Similarly, if we consider the case of two spherical resonators, as depicted in Figure 3, then we can use a bispherical coordinate system to derive explicit formulas for the capacitance coefficients [12]. Lemma 2.15. Suppose that D = D 1 ∪ D 2 consists of two spheres of radius r 1 and r 2 separated by a distance d. Define the quantities α, ξ 1 and ξ 2 as Then, it holds that In the case of larger systems of resonators we cannot hope to find such concise representations for the capacitance coefficients. However, a very useful property is that if we multiply elements in some domain B ⊂ R 3 by some factor a ∈ R, then a scaling argument can be used to see that Cap aB = aCap B . With this in mind, we can obtain explicit expressions in the case when the resonators are small compared to the distance between them. The following lemma follows from appropriate scaling arguments, where we fix the resonators and scale the distances between them [9]. Lemma 2.16. For j = 1, . . . , N , let B j be fixed, bounded subsets of R 3 with boundary in C 1,s for some 0 < s < 1. Then, consider a dilute system of N subwavelength resonators with size of order ε, given by where 0 < ε 1 and z j ∈ R 3 are fixed vectors that describe the relative position of each resonator. In the limit as ε → 0 the asymptotic behaviour of the capacitance matrix is given by Remark 2.17. It is interesting to compare and contrast the capacitance formulation to the tight-binding approximation that is commonly employed in quantum-mechanical settings [41,65]. Both these formulations construct matrix eigenvalue problems as discrete approximations to continuous differential problems. A crucial difference, however, is that the generalized capacitance matrix accounts for strong interactions between the resonators. The analogy with the tight-binding model is closer when the resonators are dilute. In this case, we see from Lemma 2.16 that C ij only depends on the i th and j th resonator and is unaffected by remaining resonators. Moreover, in the dilute regime, the eigenmodes of the system can be approximated by a linear combination of the eigenmodes of the individual resonators. This property, which is a key assumption in the tight-binding approximation, does not hold in the case of non-dilute subwavelength resonators.

Modal decompositions
The solution to the resonance problem, given in Theorem 2.7 and Corollary 2.9, can be used to understand the scattering behaviour of D. That is, we can use an expansion in terms of the resonant modes u 1 , . . . , u N (i.e. a modal decomposition) to express the scattered field when D is illuminated by some incident wave u in . We therefore consider the problem u − u in satisfies the Sommerfeld radiation condition. (2.17) Here, u is the total field while u − u in is the scattered field. We assume that the incident field satisfies ∆u in + k 2 u in = 0 in R d and ∇u in D = O(ω). The next result, from [5], shows the modal decomposition approximation of the scattered field.
Theorem 2.18. Let V be the matrix of eigenvectors of C. If ω = O( √ δ) as δ → 0 and |ω − ω i | > K √ δ for i = 1, ..., N , for some constant K > 0, then the solution to the scattering problem (2.17) can be written, uniformly for x in compact subsets of R 3 , as for coefficients a n = a n (ω) which satisfy the problem x outside of the resonators.

Higher-order approximations
The arguments used to derive the asymptotic formula in Theorem 2.7 can be continued to higher orders. For details, see [5]. For simplicity, we assume that the material parameters on each resonator are the same.
Theorem 2.20. Let d = 3. Consider a system of N subwavelength resonators in R 3 . Suppose that the material parameters are the same on each resonator, i.e. v 1 = v 2 = · · · = v N and δ 1 = δ 2 = · · · = δ N . As δ → 0, the N subwavelength resonant frequencies satisfy the asymptotic formula where λ n for n = 1, . . . , N are the eigenvalues of the generalized capacitance matrix C and τ n are given by with C being the capacitance matrix, J the N × N matrix of ones, v n the eigenvector associated to λ n and we use the norm Remark 2.21. If the material parameters v 1 , . . . , v N and δ 1 , . . . , δ N are real, then λ n and τ n from Theorem 2.20 are all non-negative real numbers. This follows from the fact that the capacitance matrix C is symmetric and positive definite. Thus, in this case the O(δ 1/2 )-term is the leading-order approximation of the real part while the imaginary part appears at O(δ).
Remark 2.22. Due to the loss of energy (e.g. to the far field), the resonant frequencies will have negative imaginary parts when the material parameters are real. In many cases it will hold that τ n = 0 for some n, meaning that the imaginary parts exhibit higher-order behaviour in δ. For example, the imaginary part of the second (dipole) frequency for a pair of identical resonators with real parameters is known to be O(δ 2 ) [17].

Two-dimensional models
Throughout Section 2 we have mainly considered the problem of a resonator array in R 3 . This was convenient for two reasons. Firstly, for small frequencies the Laplace single layer potential S 0 D approximates the Helmholtz single layer potential S ω D at leading order, in the sense that S ω D = S 0 D + O(ω) in the operator norm as ω → 0. On top of this, the fact that S 0 D is invertible in three dimensions was central to our definition of the capacitance matrix. If we consider a Helmholtz problem in two dimensions, however, we do not have either of these helpful properties. In two dimensions, S 0 D is not generally injective and the low-frequency expansion of S ω D is given by where I ∂D is the map defined as I ∂D [ϕ] = ∂D ϕ dσ for ϕ ∈ L 2 (∂D) and γ = lim n→∞ ( n k=1 1 k − log n) ≈ 0.577 . . . is the Euler-Mascheroni constant. That is, the leading-order term in the expansion of S ω D has a log ω singularity as ω → 0.
The invertibility of the Laplace single layer potential in two dimensions can be readily fixed. Let L 2 0 (∂D) be the mean-zero subspace of L 2 (∂D) defined as Then, we have that the Laplace single layer potential is well-behaved on L 2 0 (∂D) [4].
Lemma 2.23. Let d = 2. The Laplace single layer potential S 0 D is invertible from L 2 0 (∂D) onto its image. With this in hand, we can show that an invertible version of S 0 D can be defined by adding a term proportional to the integral operator I ∂D [ϕ] = ∂D ϕ dσ. For example, we have the following lemma [4].
Lemma 2.24. Let d = 2 and I ∂D be the integral map given by I ∂D [ϕ] = ∂D ϕ dσ. For any ω ∈ C \ {z ∈ C : z = iy for some y ≥ 0}, the operator S ω D : L 2 (∂D) → L 2 (∂D), defined as Notice, finally, that S ω D is nothing other than the leading-order approximation of S ω D from (2.18). This means that, up to some technical modifications, we can repeat the argument used to derive Theorem 2.7 to obtain the following analogous result.
Further to this, if for any ω and δ we define the N × N -matrix A , then the subwavelength resonant frequencies are such that the determinant of A (2) ω,δ vanishes, at leading order:

Numerical approaches
The resonant frequencies of (2.1) can be computed numerically in a variety of ways. For example, we can make use of the boundary integral formulation (2.14) to derive a discrete version of the problem. This can be achieved, for example, by discretizing the boundary using boundary elements or a multipole expansion. In order to find the resonant frequencies for a given value of δ, one needs to use a numerical root finding algorithm to find ω such that the boundary integral equation (2.14) is satisfied. The discrete version of the boundary integral operator (which depends non-linearly on ω) will need to be recomputed at each step in this iterative algorithm.
The generalized capacitance matrix can be used to obtain accurate numerical approximations with a significant reduction in computational power. Provided that δ is sufficiently small, Theorem 2.7 and Theorem 2.20 give an approximation of the resonant frequencies that is sufficiently accurate for many purposes. In Figure 4, we show the resonant frequencies of a system of ten spherical resonators computed using both the full multipole method and using the approximation with the capacitance matrix. The values derived from the eigenvalues of the capacitance matrix give a good approximation and required just 0.02 seconds of computation time, compared to the 41 seconds required for the full multipole computations, on the same computer. If greater precision than that of the capacitance matrix approximation is required, then the values derived from the capacitance matrix can be used as initial values for root finding algorithms to reduce computational time. The subwavelength resonant frequencies of a system of ten spherical resonators. We compare the values computed using the multipole expansion method to discretized the full boundary integral equation and the values computed using the capacitance matrix. The computations using the full multipole method took 41 seconds while the approximations from the capacitance matrix took just 0.02 seconds, on the same computer. Each resonator has unit radius and we use δ = 1/5000.

Periodic systems
Next, we will investigate the case when the resonators are repeated periodically as illustrated in Figure 5, often referred to as a metamaterial. Conceptually, there are three different cases depending on the dimensions of periodicity of the lattice, d l , and of the ambient space d: • d − d l = 0, which we will refer to as a crystal. In this fully periodic case (either d = d l = 2 or d = d l = 3), the metamaterial has no boundary to the surrounding space.
• d − d l = 1, which we will refer to as a screen. In this case, the structure consists of a thin sheet of resonators.
• d−d l = 2, which we will refer to as a chain. There is one example of this case, namely d = 3, d l = 1.
We assume that D, as defined in Section 2, is repeated in a periodic lattice Λ. We let P l : R d → R d l be the projection onto the first d l coordinates, and P ⊥ : We let l 1 , ..., l d l ∈ R d denote lattice vectors generating the lattice Λ, in other words such that For simplicity, we assume that P ⊥ l i = 0, which means that the lattice is aligned with the first d l coordinate axes. For a point x ∈ R d , we will throughout use the notation x = (x l , x 0 ), where x l ∈ R d l is the vector along the first d l dimensions and x 0 ∈ R d−d l . Denote by Y ⊂ R d a fundamental domain of the given lattice. Explicitly, we take The dual lattice of Λ, denoted Λ * , is generated by α 1 , ..., α d l satisfying α i · l j = 2πδ ij and P ⊥ α i = 0, for The periodically repeated i th resonator D i and the full periodic structure D are given, respectively, by v2 · · · Figure 5: A periodic array of material inclusions. Here, three material inclusions (resonators) are sketched with periodicity in one dimension. Each interior has a different wave speed v1, v2, v3 and the surrounding medium has a wave speed v. The contrast between the i th resonator and the background is given by δi, where a small value of δi describes a large contrast.
We study the equation satisfies the outgoing radiation condition as |x 0 | → ∞.
The radiation condition depends on the dimensionality (see, for instance, [43,2,16]). Note that we only impose a radiation condition for x away from the periodic structure and, in the fully-periodic case, we do not assume any radiation condition. As we shall see, since there is no radiation condition in the first d l coordinate dimensions, the spectrum σ of (3.1) is in general continuous. In order to effectively study this equation, we will use the Floquet-Bloch theory, which is outlined below.

Floquet-Bloch theory
Given a function f ∈ L 2 (R d ), the Floquet transform of f is defined as , with inverse given by (see, for instance, [16,50]) The α-quasiperiodic radiation condition depends on the dimensionality, and we refer (for example) to [16,29,31,32] for details. The spectrum σ of the original problem (3.1) is parametrised by the spectra σ(α), α ∈ Y * , of the problem (3.3), which in turn are known to consist of discrete values ω = ω α i : Definition 3.1 (Band function). The resonant frequencies ω α i , seen as functions of α, are called band functions. The collection of band functions is called the band structure.
If the spectrum σ is real, we define a band gap of D as a connected component of R \ σ, which consists of intervals in R.
As in Section 2, we will focus on the subwavelength part of the spectrum, which are the resonant frequencies ω α i which tend to 0 as δ → 0. In the case k = |α + q| for all q ∈ Λ * , we can define the quasiperiodic Green's function G α,k (x, y) as the Floquet transform of G k (x, y) in the first d l coordinate dimensions, i.e., The series in (3.7) converges uniformly for x and y in compact sets of R d , x = y, and k = |α + q| for all q ∈ Λ * . Shortly speaking, using this Green's function we can define analogous quantities and get analogous results as in the finite case. The quasiperiodic single layer potential S α,k D is then defined as The quasiperiodic single layer potential satisfies many conceptually similar properties as the "regular" single layer potential (which is a consequence of the fact that the singularity of corresponding Green's functions are the same, i.e. G α,k − G k is a smooth function of x around the origin). For example, the quasiperiodic single layer potential satisfies the jump relations on ∂D: and ) * is the quasiperiodic Neumann-Poincaré operator, given by Above, we assumed that k = |α + q| for all q ∈ Λ * . When k is small and lies in the subwavelength regime, this condition separates into two cases: • k < inf q∈Λ * |α + q|. Waves in this regime are exponentially decaying away from the structure. Such waves, which vanish in the far-field, are known as evanescent waves.
• |α| < k < inf q∈Λ * \{0} |α + q|. Waves in this regime are propagating far away from the structure, and this regime is known as the first radiation continuum. Figure 6: Example of the subwavelength band structure of a resonator array with two resonators in the unit cell. The shaded region is the first radiation continuum, defined by |α| < ω/v < inf q∈Λ * \{0} |α + q|, while the unshaded region correspond to evanescent modes. Here we see an example of band gap opening: there is an uncovered interval between the first and the second bands. Moreover, between the subwavelength bands and the higher (non-subwavelength) bands there will always be a band gap.
If we briefly assume that v i and δ i are real we can interpret the two regimes as follows. When k < inf q∈Λ * |α + q|, the problem (3.3) can be viewed as the spectral problem for a self-adjoint operator, and the resonant frequencies are real. When |α| < k < inf q∈Λ * \{0} |α + q|, due to the radiation condition this equation no longer correspond to a self-adjoint operator. Therefore, the resonators attain a small but non-zero imaginary part corresponding to the coupling with the far-field. The transition between these two regimes occurs when k = |α| for some q (known as a Rayleigh-Wood anomaly), which are the points where the spectrum becomes real and the modes become localized to the structure. In the two regimes mentioned above, we have the following integral representation (analogously to Lemma 2.1).

Evanescent-mode resonances
If we assume |α| > c > 0 for some constant c independent of ω and δ, the quasiperiodic Helmholtz single layer potential is well approximated by the corresponding Laplace single layer potential in the sense that S α,ω D = S α,0 D + O(ω 2 ) in the operator norm. Because of this, we can use analogous methods as those outlined in Section 2, and obtain similar results. It is worth pointing out that the following results hold even in the case d = 2, and for general d l .
Remark 3.7. The above theorem describes only the subwavelength part of the band structure. For small enough δ, there will be a band gap between the first N bands (which are in the subwavelength regime) and the higher bands (which are not close to 0 for small δ; see Figure 6).
Definition 3.8 (Quasiperiodic capacitance matrix). Assume α = 0. For a system of N ∈ N resonators D 1 , . . . , D N in Y we define the quasiperiodic capacitance matrix C α = (C α ij ) ∈ C N ×N to be the square matrix given by Lemma 3.9. The quasiperiodic capacitance matrix is a Hermitian matrix.
Definition 3.10 (Generalized quasiperiodic capacitance matrix). Assume α = 0. For a system of N ∈ N resonators D 1 , . . . , D N in Y we can define the generalized quasiperiodic capacitance matrix, denoted by C α = (C α ij ) ∈ C N ×N , to be the square matrix given by The next result characterizes the first N resonances of the periodic structure, and shows that they are in the subwavelength regime.
Theorem 3.11. Let d ∈ {2, 3} and 0 < d l ≤ d. Consider a system of N subwavelength resonators in Y , and assume |α| > c > 0 for some constant c independent of ω and δ. As δ → 0, the N subwavelength resonant frequencies satisfy the asymptotic formula where {λ α n : n = 1, . . . , N } are the eigenvalues of the generalized quasiperiodic capacitance matrix C α ∈ C N ×N , which satisfy λ α n = O(δ) as δ → 0. Corollary 3.13. Let d ∈ {2, 3}, 0 < d l ≤ d and assume |α| > c > 0 for some constant c independent of ω and δ. Let v α n be the eigenvector of C α associated to the eigenvalue λ α n . Then the resonant mode u α n associated to the resonant frequency ω α n is given, as δ → 0, by Figure 7 shows the resonant mode of a square crystal in two dimensions (d = d l = 2) for α close to the corner of the Brillouin zone. Here, we can observe the two-scale behaviour associated to subwavelength metamaterials: the resonant modes are oscillating on the small scale, with amplitudes which satisfy large-scale oscillations (for more details on this, we refer to [23]).

Resonances in the first radiation continuum
Here we study the regime |α| < k < inf q∈Λ * \{0} |α+q|. Since the analysis depends on the dimensionalities, we will exemplify it in the case of a metascreen, d − d l = 1.
Recall that k = ω/v. In the current regime, we must have α → 0 as ω → 0. Therefore, we assume that for some α 0 , independent of ω and such that |α 0 | < 1/v. In scattering problems, this limit corresponds to incident waves with a fixed direction of incidence (specified by α 0 ) and a frequency ω in the subwavelength regime.

Green's function and capacitance matrix formulation
In the current setting, the quasiperiodic Green's function admits the spectral representation where k 0 = k 2 − |α| 2 . The series in (3.7) converges uniformly for x in compact sets of R d , x = 0, and |α| < k < inf q∈Λ * \{0} |α + q|. In the case when k = α = 0, we define the periodic Green's function G 0,0 as When ω → 0, we then have the asymptotic expansion In particular, the Green's function has a singularity when ω → 0. In fact, this will make the analysis conceptually similar to the case in Section 2.5. We define the operator S α,k D : We then have the asymptotic expansion S ωα0,ω D = S ωα0,ω D + O(ω), with respect to the operator norm, as ω → 0. Lemma 3.14. S 0,0 D is invertible from the mean-zero space L 2 0 (∂D) onto its image.
Crucially, Lemma 3.15 shows that inverse (S ωα0,ω D ) −1 does not have the ω −1 -singularity around ω = 0. This will allow us to define capacitance coefficients in the current setting.
Definition 3.16 (Periodic capacitance matrix). For α 0 with |α 0 | < 1, and for a system of N ∈ N resonators D 1 , . . . , D N in Y we can define the periodic capacitance matrix C 0 = (C 0 ij ) ∈ R N ×N to be the square matrix given by Similarly to Lemma 2.11, the periodic capacitance matrix can alternatively be written Because of this formula, C 0 share many of the properties of the capacitance matrices in previous settings.
Lemma 3.17. The periodic capacitance matrix C 0 is a real, symmetric, positive semi-definite matrix with one vanishing eigenvalue. Moreover, C 0 is independent of α 0 .
The name periodic comes from the fact that an α-quasiperiodic function is, in the case α = 0, a periodic function. Since α = ωα 0 and ω is small, we are working close to the periodic case. In fact, since C 0 is independent of α 0 , all relations involving C 0 are equal to the periodic case α 0 = 0. Definition 3.18 (Generalized periodic capacitance matrix). For a system of N ∈ N resonators D 1 , . . . , D N in Y we can define the generalized periodic capacitance matrix, denoted by C 0 = (C 0 ij ) ∈ C N ×N , to be the square matrix given by Since S ωα0,ω D has a ω −1 -singularity as ω → 0, we cannot apply the Gohberg-Sigal theory to study the characteristic value perturbation of A α (ω, δ) defined in Lemma 3.4. Instead, we rephrase the system in the following result [11].
Lemma 3.19. Let k = ω/v and assume that α = ωα 0 for some α 0 independent of ω and δ such that Now we can apply the functional analytic approach, outlined in Section 1.3, to derive the following result on the subwavelength resonant frequencies in the current setting [11,8].

Higher-order approximations
In the case δ i v 2 i ∈ R for all i, the leading-order parts of the resonances, given in Theorem 3.20, are real. It is also interesting to compute corresponding imaginary parts, which will specify the bandwidth of the resonant behaviour. We have the asymptotic expansion of S ωα0,ω . Then we can define the "higher-order" matrix C 1,α0 as and corresponding generalized matrix C 1,α0 as Therefore, we have the following result, which gives the resonances ω α i to a higher order [11]. Theorem 3.21. Let d − d l = 1, k = ω/v and assume that α = ωα 0 for some α 0 independent of ω and δ such that |α 0 | < 1/v. As δ → 0, the N subwavelength resonant frequencies satisfy ω α n = ω α n + O(δ 3/2 ) where ω α n , for n = 1, ..., N , are the roots ω = ω α n of the equation While the capacitance coefficients can be thought of as total charge (or "mass"), c i is the centre of mass (up to rescaling). Briefly put, we can compute C 1,α0 in terms of these coefficients, which allows us to compute explicit expressions of ω α n .

Modal decomposition
The "higher-order" term S α0 −1 in the expansion of the single layer potential enters the expression for the resonant modes in this case. From the arguments used to derive Theorem 3.21 we have the following result on the resonant modes [11].
Theorem 3.22. Let d − d l = 1, k = ω/v and assume that α = ωα 0 for some α 0 independent of ω and δ such that |α 0 | < 1/v. Let v 0 n be the eigenvector of C 0 associated to the eigenvalue λ 0 n . Then the resonant mode u α n associated to the resonant frequency ω α n is given, as δ → 0, by where S α,k D : R d → C N is the vector-valued function given by . With these resonant modes at hand, we can easily solve the associated scattering problem u α − u in satisfies α-quasiperiodic radiation condition as |x 0 | → ∞. Here, we assume that u in is a plane wave with some wave vector k; u in (x) = e ik·x . Furthermore, α is now specified in terms of k as α = P l k, where (as before) P l is the projection to the first d l coordinates.
Then, for small ω we assume that k = ωw for some fixed w.
Theorem 3.23. Let d − d l = 1, and assume that u in = e ik·x where k = ωw. Also, assume that |ω − ω 0 i | > K √ δ for i = 1, ..., N , for some constant K > 0. Then, as δ → 0, the solution to the scattering problem (3.12) can be written, uniformly for x in compact subsets of R d , as

13)
where V is the matrix of eigenvectors of C 0 and the coefficients a n = a n (ω) satisfy the problem

Connection to quasiperiodic capacitance matrix
In the previous sections, we have used the notation C 0 for the periodic capacitance matrix, and C α for the quasiperiodic capacitance matrix, defined for α ∈ Y * \ {0}. This choice of suggestive notation is deliberate, and in this section we will show that these capacitance matrices can be combined into a continuous function of α ∈ Y * , in the sense that C 0 = lim α→0 C α . For α ∈ Y * , α = 0, the Green's function G α,0 satisfies so that for small α we have In this setting we define the operator S α,0 D : L 2 (∂D) → H 1 (∂D) as The structure of this operator is entirely analogous to S α,k D defined in Section 3.3.1 (where, conceptually, the limit ω → 0 now corresponds to |α| → 0). Therefore, we can apply the same method to obtain the following result.
Lemma 3.24. The periodic capacitance matrix C 0 and the quasiperiodic capacitance matrix C α , for

Applications of the generalized capacitance matrix
The generalized capacitance matrix can be applied to explain a variety of interesting physical phenomena. We particularly want to study extraordinary macroscopic properties (such as exotic effective parameters, rainbow trapping or unidirectional scattering) as well as robust localization at subwavelength scales.

Double-negative materials
The field of metamaterials more or less began with the realisation that micro-structured media could be designed to have effectively negative material parameters. This behaviour can also be replicated in the setting studied in this work, as reported in [17]. Consider a large number N of identical resonator pairs, given by where D = D 1 ∪ D 2 is some fixed pair of resonators, 0 < s 1 is some characteristic size, z N j ∈ R 3 are the positions of the resonator pairs and R d N j are rotations in R 3 which orient the pair in the direction of the unit vector d N j . We will assume that sN = Λ for some fixed Λ > 0 and that there is some bounded domain Ω such that {z N j : 1 ≤ j ≤ N } ⊂ Ω for any N ≥ 1. We also want that the resonators are regularly distributed in the sense that there exists some number ν such that Additionally, we want that there exists some positive function V ∈ C 1 (Ω) and a matrix-valued function B ∈ C 1 (Ω) such that there are constants C 1 and C 2 which satisfy 3) for all f ∈ C 0,α (Ω) with 0 < α ≤ 1 and all N ≥ 1, and similarly  We must make some additional assumptions to achieve the desired double-negative behaviour. We assume that the material parameters are such that v 1 = v 2 ∈ R and δ 1 = δ 2 = µ 2 s 2 for some number µ > 0. This means that the generalized capacitance coefficients associated to a single resonator pair do not depend on s. As such, the leading-order term in the expansion of the subwavelength resonant frequencies, as s → 0, is fixed. Additionally, we will assume that the pair of resonators is symmetric in the sense that for ψ 1 and ψ 2 as defined in Lemma 2.2. In this case, we also have that the subwavelength resonant frequencies for the resonator pair are given, as s → 0, by for real numbers τ 1 , η 1 and η 2 .
We want the incident frequency to be close to the second resonant frequency, which corresponds to the dipole resonant mode. In particular, we assume that there is some a < µ 3 η 1 such that ω = λ 2 + as 2 . (4.8) Under these assumptions, we can derive an effective medium theory for the system as N → ∞. Let u N be the field scattered by an array of N resonator pairs of the form (4.1). Suppose there is some macroscopic field u ∈ C 1,α (Ω) which is such that u N converges to u in C 1,α (Ω) as N → ∞. Then, we can show that u must satisfy the equation [17] ∇ where for positive constants g 0 and g 1 given by Equation (4.9) shows the desired double-negative behaviour. Since we assumed that a < µ 3 η 1 , it holds that g 1 > 0. Thus, if B(x) is a positive definite matrix for all x ∈ Ω and Λ is chosen to be sufficiently large, then M 1 (x) will be a negative definite matrix for all x ∈ Ω. Likewise, since g 0 > 0 and V > 0 by assumption, we can choose Λ to be sufficiently large that M 2 (x) < 0 for all x ∈ Ω. We refer to [17] for the details.

Frequency separation with graded arrays
By introducing a gradient to the properties of a metamaterial, it is possible to cause different frequencies to be localized at different positions in the structure. This frequency separation is useful for a variety of applications since it means an incoming signal can be separated into its different frequency components. This phenomenon is often known as rainbow trapping and has been observed in a variety of electromagnetic [62,46], acoustic [4,69,47] and other metamaterials [30]. Approximating the system using the generalized capacitance matrix, we can study the extent to which a system of coupled subwavelength resonators is able to perform frequency separation in a controllable manner.
This behaviour is very similar to the action of the cochlea. Devices based on these principles have been used to design biomimemtic hearing devices [6,59,49]. For example, a graded resonator array is shown in Figure 8 to replicate the frequency separation of the cochlea. These structures are useful for building artificial hearing approaches as well as, conversely, learning about the function of human hearing. For example, understanding the details of cochlear amplification is a significant open question and is obscured by the challenges in experimenting on living biological organisms. By designing analogue artificial devices, we are able to test theories and reveal crucial insight into the fundamantal mechanisms that underpin human hearing [6,59].

Sensitivity enhancement using high-order exceptional points
We want to design "enhanced" sensors that are strongly influenced by small perturbations in their immediate surroundings. For example, we want to be able to sense the presence of a small object, as depicted in Figure 9a. This small object might be a virus, for example [64]. Typically, such an occurrence would cause a shift in the resonant frequencies that is proportional to the strength of the perturbation (i.e. proportional to the size of the small object). The idea here, however, is to design an array for which this shift is enhanced, particularly for small perturbations [66,45]. We achieve this by adding sources of energy gain and loss to the system, represented by material parameters with non-zero imaginary parts. This means that Lemma 2.13 no longer holds and we can fine-tune the material parameters in order to create systems with coincident eigenvectors. The fundamental result guiding this work is the following lemma [9]. Lemma 4.2. Let m ∈ L ∞ (R 3 ) and N ∈ N. Suppose that a system has an N th -order singularity in the sense that the Green's function G k m , which is defined as the solution to ∆ (a) A graded array of resonators. High frequencies will give a maximum response towards the left of the array while lower frequencies will be detected further to the right.  PT -symmetric array of active resonators with an exceptional point small object (a) An active system of resonators can be designed such that it supports an exceptional point, meaning that it experiences an enhanced response to small perturbations such as the presence of a small particle.  for k in a neighbourhood of some k * ∈ C, where G k is the free space Helmholtz Green's function, {ϕ 1 , . . . , ϕ N } ⊂ H 1 loc (R 3 ) are generalized eigenfunctions associated to k * and the remainder R is a holomorphic function of k that is smooth as a function of x and y. If a small material inclusion Ω is introduced to the system, then the new system has a resonant frequency k Ω with the asymptotic behaviour with η z is a constant that depends on the position and material parameters of the small particle Ω, the background material parameters m and the generalized eigenfunction ϕ N .
One way to create N th -order singularities, as are required by Lemma 4.2, is to design structures with higher-order resonant singularities. That is, structures with exceptional points, where eigenvalues and eigenvectors coincide. In this setting, we will search for asymptotic exceptional points. Remark 4.4. The restriction to considering asymptotic exceptional points in this work is not a weakness of the analytic approach but represents the behaviour of the system. The radiation condition means that the symmetry we impose on the resonators is not extended to the far field, meaning that we won't have exact degeneracy at the exceptional points.
Exceptional points are a consequence of balanced symmetries in the system, which cause the eigenvectors to align. So that the system already has some underlying symmetry, exceptional points are often sought in structures with parity-time symmetry. We will assume that the problem is parity-time symmetric in the sense that each resonator D i can be uniquely associated to another resonator D j (possibly with i = j) such that where the parity operator P : R 3 → R 3 and the time-reversal operator T : C → C are given, respectively, by P(x) = −x and T (z) = z.
For a system of two subwavelength resonators, we can find the eigenvalues and eigenvectors of the generalized capacitance matrix explicitly. Using these formulas, we can show that there exists an asymptotic exceptional point for certain parameter values [11].
The approximate nature of the asymptotic exceptional point predicted in Theorem 4.5 is demonstrated by Figure 10, where the subwavelength resonant frequencies of the full differential system are simulated directly using the multipole expansion method. We can see that there is a critical value of the gain and loss (the imaginary parts of the material parameters) such that the eigenvalues coincide at leading order.  Below this critical value the leading-order parts of the resonant frequencies are real and they form a purely imaginary conjugate pair above this value. In order to produce higher-order subwavelength exceptional points, we need to study larger systems of resonators. So that the capacitance matrix is easier to work with in this case, we will make an additional assumption that the resonators are relatively far apart, meaning that we can use the dilute approximation given in Lemma 2.16. This means we can efficiently analyse large systems of resonators, in terms of this leading-order approximation of the capacitance matrix [9]. For a system of three subwavelength resonators, we can show that there is one third-order subwavelength asymptotic exceptional point. This is described by Theorem 4.6 and the subwavelength resonant frequencies are depicted in Figure 11. Theorem 4.6. A PT -symmetric system D of three dilute resonators has an asymptotic exceptional point of order 3 with respect to ε and δ at the resonant frequency ω * , which is given as ε, δ → 0 by where c 1 is the real root of the polynomial c 3 1 + 27 4 c 1 − 27 8 = 0 (i.e. c 1 ≈ 0.483...). We can continue this process to study higher-order exceptional points in larger systems. We quickly find that as the dimension grows the number of exceptional point similarly grows. In a system of four subwavelength resonators we find that there are four asymptotic exceptional points. The imaginary parts of the material parameters on each resonator at each exceptional point are depicted in Figure 12. We can see that each exceptional point corresponds to one of the four different combinations of relative magnitude and sign that is possible under the assumption of PT symmetry.
The symmetry exhibited by the fourth-order asymptotic exceptional points shown in Figure 12 can also be seen in higher-order exceptional points in larger structures. For example, in Figure 13, we take the fourth-order exceptional points from Figures 12a and 12d and find asymptotic exceptional points of order 8 and 14 in dilute resonator arrays of the corresponding size. Exceptional points with the same qualitative distribution as in Figures 13d-13f were previously observed in Hamiltonian systems in [68]. This analysis demonstrates the value of the generalized capacitance matrix, particularly under an assumption of diluteness, which gives a concise yet rigorous approximation of the behaviour of the differential system.

Subwavelength guiding of waves
The study of point defects in band gap materials has had immense impact on technological applications. The most notable examples are doped semiconductors, where conducting modes are induced by impurities   Figure 11: A PT -symmetric system of three subwavelength resonators supports an asymptotic exceptional point, where the eigenvalues (and corresponding eigenvectors) of the dilute capacitance matrix coincide. We can compare the resonant frequencies of the full differential problem in (a), computed using the multipole expansion method, and the approximate frequencies in (b), computed using the dilute approximation of the generalized capacitance matrix. In both cases, the relative size of the imaginary parts are fixed and they are rescaled by some τ , which is such that the asymptotic exceptional point occurs at τ = 1.  12: A system of four PT -symmetric resonators supports four asymptotic exceptional points, each with its own symmetry. Here, we plot the leading order coefficients of the imaginary parts of the material parameters (the gain or loss) at each of the four asymptotic exceptional points, using the notation b = Im(v 2 1 δ1) and d = Im(v 2 2 δ2). in the semiconductor material. As we have seen, band gaps can be found in any type of wave-propagation systems, and defects in a subwavelength band gap material can be used to enable trapped or guided waves on very small length-scales.
To fix the setting, we consider the equation in the fully periodic case d = d l = 2. We study the equation (3.1), where D consists of a single, circular resonator inside Y . However, instead of considering the periodic crystal D, we detune the size of some resonators, thereby creating a defect. We will consider the two defects illustrated in Figure 14, where either a single resonator or a line of resonators are detuned. We choose a square lattice with unit cell Y and let D be a circle of radius R and D d a circle of radius R + ε for some −R < ε < 1 − R. Then we define the defect crystals We will always consider defect structures in relation to a corresponding unperturbed, periodic structure, in this case D.

Fictitious source method for point defects
We wish to solve (3.1) where the periodic crystal D is replaced by D pt or D ln . Firstly, we observe that D ln is periodic in the x 1 -direction, so corresponding problem can be reduced to the quasiperiodic problem in the strip Y strip = − 1 2 , 1 2 × R by Floquet-Bloch theory. This makes the analysis of the two problems (D pt on R 2 or D ln on Y strip ) conceptually very similar. We will outline the method of fictitious sources in the case of a point defect. We are studying the problem where k r = ω/v r and v r is the wave speed inside the resonators. The idea is to replace the defected resonator with an unperturbed resonator, along with fictitious sources which are designed to make the new problem equivalent to the original. More precisely, we consider the problem where f, g are the source terms and χ ∂D is, as usual, the characteristic function of ∂D. Note that (4.12) is posed on the periodic geometry D, and that the non-zero sources are present only on the boundary of the central resonator D. We seek solutions u of (4.11) and u of (4.12) for ω inside the band gap of the periodic problem. Inside the unit cell Y , we can represent the solution u as in D, (4.13) where H satisfies the homogeneous equation ∆H + k 2 H = 0 in Y . Using this approach, we can explicitly compute f, g in terms of ψ, φ so that u coincides with u in Y \ (D ∩ D d ).
Lemma 4.8. The density pair (ψ, φ) and the effective sources (f, g) satisfy the following relation where A ε is defined as Here, A and A D d are the operators defined in Lemma 2.1 for the domains D and D d , respectively, and the operators P 1 : L 2 (∂D) 2 → L 2 (∂D d ) 2 and P 2 : L 2 (∂D) 2 → L 2 (∂D d ) 2 are defined by As in previous settings, we can use an integral equation formulation of the problem. In the defect crystal setting, however, the formulation is slightly different than previously [15].
Lemma 4.9. The subwavelength band-gap frequencies of (4.11) are precisely the characteristic values ω = ω ε (δ) of the operator inside the band gap of D, such that ω ε → 0 as δ → 0. Here, A α is the operator defined in Lemma 3.4, which is invertible for ω inside the band gap.
We let ω * = ω * (δ) be the maximum of the first band: The following theorem demonstrates the existence of band-gap frequencies in the point defect crystal [15].
Theorem 4.10. Assume that δ and ε are small enough and the pair (R, ε) satisfies one of the two assumptions: (i) R small enough and ε < 0 (dilute regime); (ii) R close enough to 1/2 and ε > 0 (non-dilute regime).
Then there exists a subwavelength band-gap frequency ω ε of (4.11). In both cases we have the asymptotic expansion when ε and δ go to zero, for some constant µ > 0.
Remark 4.11. An explicit expression for µ is given in [15].
Remark 4.12. In the reverse cases, i.e. if R is close to 1/2 and ε < 0 or if R is small and ε > 0, there will be no band-gap frequencies. In these cases, the resonance frequency of the defected resonator is shifted downwards, and will therefore lie outside of the band gap.

Guided waves in line defects
We let D ln be the line defect crystal and study the problem In this case, we have a very similar characterization of the band-gap frequencies to the one given in Lemma 4.9.  Lemma 4.13. The subwavelength band-gap frequencies of (4.11) are precisely the characteristic values ω = ω ε (δ, α 1 ) of the operator inside the band gap of D, such that ω ε → 0 as δ → 0.
We now let ω α1, * be the maximum of the first band at α = (α 1 , α 2 ): Then we have the following result.
Theorem 4.14. Assume that δ and ε are small enough and the pair (R, ε) satisfies one of the two assumptions: (i) R small enough and ε < 0 (dilute regime); (ii) R close enough to 1/2 and ε > 0 (non-dilute regime).
In order for the whole defect band to lie in the subwavelength band gap, we need a sufficiently large ε. Since Theorem 4.14 is based on asymptotically expanding M for small ε, a different analysis is needed to handle this case. The following result is based on asymptotics in the dilute regime, and is valid even for ε with large magnitude [21].
Proposition 4.16. For R and δ small enough, there exists an ε 0 > 0 such that for any ε ∈ (−R, −ε 0 ) we have In order for the line defect crystal to be useful as a waveguide, we need the localized modes to propagate along the defect line. In other words, we must exclude the case of bound modes, which are modes that are localized along the direction of the line. If there is such a mode u, corresponding to a frequency ω, we can apply the Floquet transform so that u α (x) := F[u](x, α) solves (4.11) for any α. Corresponding band function attains the same value ω for any α, so we conclude that bound modes are associated to flat band functions. The next result from [21] shows that the defect modes in our case are not bound along the defect line.  Here ε < −ε0, which means that the whole defect band will be above the bulk band. Moreover, we see that the defect band is not flat, and corresponding modes are propagating along the line defect. The defect band was computed using the asymptotic formula (4.20) (red dashed) and by discretizing the operator M ε,α 1 (ω, δ) (defined in (4.18)) using the multipole method (red circles).
The integral representation for the localized and guided waves (Lemma 4.9 and Lemma 4.13) can be discretized using the multipole method, which provides an efficient method to compute the band-gap frequency and corresponding localized mode (see Figure 15 and Figure 16).

Robust guiding at subwavelength scales
There is a fundamental restriction of the practical applicability of the localized modes studied in Section 4.4. Taking the point defect as example, we have from Theorem 4.10 that the band-gap frequency is exponentially close to the edge of the bulk bands (see Figure 15). When fabricating such structures, there is a large risk that small imperfections will cause the band-gap frequency to be lost inside the bulk bands, so that the desirable localization property disappears.
For structures to have physically achievable localization properties, such properties must be robust against imperfections of their design. To accomplish this, we take inspiration from the field of topological insulators. These are materials which can conduct current along the edges of the material, while are insulating in the bulk. The striking property of topological insulators is that these conducting edge modes originate from the structure of the bulk, rather than the edge itself. Moreover, they are localized to the edges and exhibit a remarkable robustness against perturbations of the system.
We study a high-contrast resonator analogue of the Su-Schrieffer-Heeger (SSH) model [61]. We consider the case of a chain of resonators, corresponding to d = 3 and d l = 1, where the unit cell is given by Moreover, we consider the case of a dimer of resonators D = D 1 ∪ D 2 . We need two assumptions of symmetry for the analysis that follows. The first assumption is that each individual resonator is symmetric in the sense that there exists some x 1 ∈ R such that where R 1 and R 2 are the reflections in the planes p 1 = {−x 1 } × R 2 and p 2 = {x 1 } × R 2 , respectively. We also assume that the dimer is symmetric in the sense that where, as before P(x) = −x. We denote the resonator separation within the cell as l, and between the cells as l , i.e. l = 2x 1 , l = L − l, see Figure 17. For simplicity, we choose v = v 1 = v 2 = 1, and study effects originating from the geometry of the structure. As we shall see, a topological phase transition occurs when l changes across the symmetry point l = L/2 (corresponding to l = l ). We begin by considering the periodic equation satisfies the α-quasiperiodic outgoing radiation condition as x 2 2 + x 2 3 → ∞. (4.23) To enable explicit computations of asymptotic expansions and topological properties, we will assume that the resonators are dilute, in the sense that they occupy a small volume compared to the surrounding medium. As in Section 4.3, we assume that the resonators can be obtained by rescaling fixed domains B 1 , B 2 as follows: for some small parameter ε > 0. The first result shows that, in addition to the band gap above the second band, there is also a band gap between the first two bands in the case l = L/2 [7].
Theorem 4.18. In the dilute regime and with δ sufficiently small, there exists a subwavelength band gap between the first two bands if l = L/2, i.e.
for ε and δ small enough. Figure 17: Example of the array, drawn to illustrate the symmetry assumptions.
In the case l = l = L/2, the first and the second bands will have a degeneracy at α = π/L. Moreover, the two bands will intersect linearly in a so-called Dirac cone. By altering l, this Dirac cone will open a band gap. As we shall see, the nature of this band gap is fundamentally different in the two cases l < l and l > l .

Topological indices and band inversion
We will now define an index which quantifies the topological properties of the Bloch eigenbundle.
Definition 4.19 (Zak phase). For a non-degenerate band ω α j , we let u α j be a family of normalized eigenmodes which depends continuously on α. We then define the Zak phase ϕ zak j as ϕ zak where ·, · denotes the L 2 (D)-inner product.
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 the current setting, a monopole mode is a mode with even parity, while a dipole mode is a mode with odd parity. For the two-resonator chain studied here, the first mode will always be of monopole nature at the origin α = 0 (which is a consequence of the fact that the eigenvectors of the periodic capacitance matrix C 0 are, respectively, ( 1 1 ) and 1 −1 ). At α = π/L, the modes will also be of monopole/dipole nature, but which of these corresponds to the first or second mode depends on the geometry of the structure [10].
Theorem 4.20. We assume that D is in the dilute regime specified by (4.24). Then the Zak phase satisfies if l > l , for ε and δ small enough.
The Zak phase predicts band inversion in the case l > l , which is the statement of the next result.
Proposition 4.21. For ε and δ small enough, the crystal has undergone band inversion in the case l > l , but not if l < l . In other words, the eigenfunctions associated with the first and second bands at α = π/L satisfy u π/L 1 , when l > l .
As δ → 0, these eigenfunctions are given, respectively, by Here, S α,ω D (x) is the function defined in Corollary 3.13.

Robustness of edge modes
We now study a finite chain of resonators which supports topologically protected edge modes. Specifically, we assume that D has the form 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 l and l that are swapped at the middle resonator.
An example of such a configuration is depicted in Figure 18. In this figure, it is shown how to associate different Zak phases on either side of the central resonator (which constitutes the "edge"). Based on the principle of bulk-boundary correspondence, we thereby expect robust localized modes around this edge. We model wave propagation in the crystal D by the Helmholtz problem (4.26) Figure 19 shows a comparison between the band-gap frequencies of a chain with a topological defect (as in Figure 18), compared against a chain with a point defect (analogously as in Section 4.4). Not only is the topological band-gap frequency further away from the edges of the band gap, but it also exhibits a lower variance when random Gaussian errors is imposed on the resonator locations. For large error, we see that the band-gap frequency of the topological defect chain (Figure 19a) is much more robust than the band-gap frequency of the point-defect chain. In Remark 2.17 we observed that, under the dilute assumption, the capacitance formulation is analogous to the tight-binding model commonly utilized in studies of quantum-mechanical systems. The tight-binding model is often coupled with a nearest-neighbour approximation, whereby long-range interactions are neglected. Figure 20b shows the bulk and band-gap frequencies computed using such nearest-neighbour approximation. Compared to Figure 20a, where all interactions are taken into account, we observe that the nearest-neighbour approximation is not accurate. This discrepancy shows a fundamental difference between topological edge modes in the setting of classical waves and quantummechanical waves. Moreover, we see that the band-gap frequency in Figure 20b is unaffected by the error. This can be seen as a consequence of chiral symmetry, which is not present without the nearest-neighbour approximation.

Edge modes in a dislocated chain
In Section 4.5.1 we saw that the structure under consideration has different topological properties in the different cases l > l and l < l . These two regimes only differ by a choice of unit cell: if we shift the unit cell Y by L/2 the regimes will be swapped. Correspondingly, edges between the different topological regimes are created by introducing a shift between the unit cells on either side of the edge. This is exemplified in Figure 18. Another example of this principle is to translate, or dislocate, half the chain by some distance d as illustrated in Figure 21. In this case, we consider a chain of resonators given by      where D is the resonator dimer defined as in the beginning of this section.
Qualitatively, it is straightforward to understand how the dislocation will affect the existence and behaviour of band-gap frequencies. As d increases from 0, we will detune the dimer coupling, which is responsible for the band gap, and band-gap frequencies will therefore appear from each edge of the band gap. On the other hand, when d is very large the two half-space chains will decouple, and the bulk-boundary correspondence suggests that each of these half-chains will support a single band-gap frequency. We therefore expect only a single band-gap frequency in this limit. As d varies between 0 and ∞ these two frequencies hybridize, and will together cover the whole band gap. This is sketched in Figure 22.
We now set out to prove that the picture outlined above is correct. Similarly to Section 4.4, we will apply a fictitious source method to model the dislocated structure D disloc in terms of the original structure D along with fictitious sources on the dislocated resonators. The analysis is separated into three parts, depending on the value of d: In this case, we can use asymptotic expansions in terms of d to prove that there is a band-gap frequency emerging from each edge of the band gap.
• d = mL for some m > 0. This particular dislocation is equivalent to removing m dimers from D.
This observation simplifies the problem and allows explicit computations of the band-gap frequencies in terms of the eigenvalue problem of certain Toeplitz matrices. With these ideas at hand, the following two theorems were proved in [7]. The first result, valid for small d, shows the emergence of a band-gap frequency from each edge of the band gap.
Theorem 4.23. Assume that the resonators are in the dilute regime specified by (4.24) and that l 0 > 1/2. Then, for small enough δ and ε, there exists some d 0 = O(ε) such that there are two band-gap frequencies ω 1 (d) and ω 2 (d) for all d ∈ [d 0 , ∞), both of which converge to the same value ω ∞ as d → ∞.
In particular, Theorem 4.23 states that the band-gap frequencies will cover an interval I := [ω 1 (d 0 ), ω 2 (d 0 )] inside the band gap, and therefore allows us to fine-tune the system to achieve optimal robustness. In Figure 23 we see numerically computed eigenmodes, clearly showing how the edge modes of the "half system" hybridize to create the two band-gap frequencies. Here, we simulate an array with parameters L = 9, l = 6, d = 30 and δ = 1/7000. Right: For comparison, the edge mode of the corresponding 'half system' is shown, which can be thought of as the d = ∞ case. In both cases, the eigenmodes are shown directly above the corresponding system of resonators.

Non-Hermitian band inversion and edge modes
In the previous sections, we assumed that the wave speeds satisfy v = v 1 = v 2 = 1, and only considered topological phenomena originating from the geometry of the resonator chain. In this section we do the opposite, namely we study topologically protected modes in structures where the geometry is periodic, and where topological edges are introduced in terms of the resonator wave speeds. Most importantly, we allow the wave speeds to be complex, giving a non-Hermitian capacitance formulation similarly as in Section 4.3. Details of this analysis are found in [18].
We let D be the dimer D = D 1 ∪ D 2 as defined previously in Section 4.5.1. For i = 1, 2 and m ∈ Z, we introduce the notation D m i := D i + (mL, 0, 0), Here, v m i denotes the wave speed inside D m i . We will assume that the contrast parameter is given by δ ∈ R for all resonators. We then consider satisfies the outgoing radiation condition as x 2 2 + x 2 3 → ∞. (4.27) In the general setting, this equation is not periodic and cannot be approached using Floquet-Bloch theory.
For prescribed values v 1 , v 2 ∈ C we will consider two different cases, depending on v m i : We emphasize that complex values of v i correspond to non-Hermitian structures. In the case v 1 = v 2 , the periodic structure is PT -symmetric.
We begin by studying the periodic structure. In Section 4.5.1 we described the concept of band inversion in the Hermitian case. In the present case, since the quasiperiodic generalized capacitance matrix C α is non-Hermitian, the band inversion phenomenon is formulated in a slightly different manner. For j = 1, 2, we let u j and w j denote a bi-orthogonal system of eigenvectors of C α . In other words, u j and w j are eigenvectors of C α and (C α ) * , respectively, and satisfy w i · u j = δ ij . We let the right Bloch eigenmodes u α j be the modes of (4.27), and the left Bloch eigenmodes w α j be the modes of (4.27) with wave speeds given by v 1 and v 2 , respectively. From Corollary 3.13, we then have Definition 4.24 (non-Hermitian Zak phase). For a non-degenerate band ω α j of the periodic structure (4.28), we let u α j and w α j be a family of normalized right, respectively left, eigenmodes which depend continuously on α. We define the (non-Hermitian) Zak phase ϕ zak In the case v 1 = v 2 ∈ R, this definition coincides with Definition 4.19, and we therefore choose the same notation for these two definitions. In the sequel, we will occasionally write ϕ zak j (v 1 , v 2 ) to denote the Zak phase defined with wave speed v 1 inside D 1 and v 2 inside D 2 .
In the Hermitian case, a non-zero Zak phase is equivalent to an inverted band structure. The fact that the Hermitian Zak phase is quantized originates from the fact that the eigenmodes are purely monopole and dipole modes at α = 0 and α = π/L. Unlike the Hermitian case, the non-Hermitian Zak phase is not quantized. A non-integer value of the Zak phase can be attributed to a "partial" band inversion, where the eigenmodes are expressed as (complex) linear combinations of monopole and dipole modes, which swap as α traverses Y * .
For the next result, we will assume that the Hermitian counterpart of the structure is topologically trivial. More precisely, we assume Then the following result from [18] holds.
Proposition 4.25. Assume that the chain is periodic and topologically trivial, in other words that it satisfies (4.28) and (4.30). Then we have Remark 4.26. Proposition 4.25 provides intuition on how to create structures supporting edge modes. Proposition 4.25 shows that distinct Zak phases can, in general, be achieved by swapping v 1 and v 2 while keeping the geometry fixed. Therefore, the defect specified in (4.29) results in Zak phases with opposite signs on the different sides of the edge. The reason we assume (4.30) is to emphasize that distinct Zak phases can originate as a pure non-Hermitian effect, which disappears in the Hermitian limit as Im(v 1 ), Im(v 2 ) → 0.
We now turn to the analysis of the defect problem, specified in (4.29) and illustrated in Figure 24. Since the geometry of this defect chain is periodic, we can utilize a Floquet-Bloch approach to derive a capacitance matrix characterization of any localized mode. We define the matrices Figure 24: Illustration of the edge. The special case v1 = v2 corresponds to local PT symmetry. Legend: material parameter v1, material parameter v2.
for some parameter b ∈ C, which can be interpreted as the decay of the localized mode between two consecutive resonators. The corresponding capacitance matrix, whose eigenvalue problem provides a discrete approximation to the localized mode and corresponding frequency, is given by The eigenvalues of this matrix depend on α ∈ Y * . In order for localized modes to exist, there must be an with λ 1 = C π/L 11 + C π/L 12 and λ 2 = 2C 0 11 . Depending on the values of v 1 and v 2 , we either have |b + | = |b − | = 1 (in which case there are no localized modes) or that |b ± | < 1 while |b ∓ | > 1. Based on these ideas, we can prove the following result [18].
Theorem 4.27. Suppose that the array of resonators has a defect in the material parameters specified by (4.29). Then, for small δ, l 2 −1 (unbroken PT -symmetry), the structure does not support simple localized modes in the subwavelength regime.
(broken PT -symmetry) or if v 1 = v 2 (no PT -symmetry), the frequency ω of a simple localized mode in the subwavelength regime must satisfy Here, b 0 is the value of b specified by (4.31) satisfying |b 0 | < 1.
Remark 4.28. Theorem 4.27 provides a characterization of the possible frequency ω and decay length b of simple localized modes in the subwavelength regime. In order to prove that such mode indeed exists, we would have to prove that one eigenvalue µ α j (b 0 ) is indeed constant in α. Analytically, this is obscured by the fact that the capacitance coefficients have a complicated dependency on α. Numerically, however, the eigenvalues µ are straightforward to compute, providing convincing evidence for such modes to exist. Moreover, the localized modes can easily be computed in large finite chains (see Figure 25), demonstrating excellent agreement between the numerical and theoretically predicted values.
In both the Hermitian and non-Hermitian cases, there are degeneracies associated to symmetric structures. Breaking the symmetry can open the degeneracy into topologically distinct band gaps as illustrated in Figure 26. In the Hermitian case, the degeneracy is a linear intersection known as a Dirac cone. In the non-Hermitian, PT -symmetric case, there can be exceptional point degeneracies which open when the PT -symmetry is detuned. In both the Hermitian and the non-Hermitian case, edge modes are created when the two topologically different phases are joined along an interface.

Bound states in the continuum and Fano resonances
Sometimes, localized modes can exist in periodic structures without a defect. This typically happens when the structure has certain symmetries, resulting in resonant modes in the radiation continuum whose far-field radiation vanishes. Such states are known as bound states in the continuum.
Definition 4.29. A resonant mode u α n of (3.3) is said to be a bound state in the continuum if the corresponding resonant frequency ω α n is real, satisfies |α| < ω α n /v and the mode satisfies, As we shall see, there are two symmetry conditions required to achieve bound states in the continuum in the subwavelength resonator arrays. The first condition is a symmetry condition of the structure D, while the second condition is that α = 0, which corresponds to modes that radiate perpendicularly to the structure. We then have the following result from [8].
Theorem 4.30. Assume that d − d l = 1 and that D = D 1 ∪ D 2 ⊂ Y . We assume that D satisfies the symmetry conditions and δ 1 v 2 2 = δ 2 v 2 2 ∈ R, where P, P 0 : R d → R d are the parity operators P(x) = −x and P 0 (x l , x 0 ) = (x l , −x 0 ).
Moreover, we assume that α = 0. Then, for small enough δ, the second resonant mode u 0 2 is a bound state in the continuum.
Interestingly, we can never have the first resonant mode u 0 1 as a bound state in the continuum. This is due to the fact that the first mode corresponds to the broad, low-frequency response of the screen which is not due to the local resonance of the structure.
An interesting problem is now to describe the behaviour when the symmetry conditions are no longer satisfied. In particular, when P 0 D i = D i , we can use Theorem 3.21 to conclude that ω α 2 has a small but non-zero imaginary part. As we will see next, this causes an interesting Fano-type transmission scattering.
Given a unit vector w = (w l , w 0 ) ∈ R d with w 0 > 0 we define the wave vectors We also let α = ω v w l . We now assume that and seek the behaviour of the solution u of (3.12). In the first radiation continuum, the scattered field u − u in consists of a single propagating mode as |x 0 | → ∞. We will write f ∼ g to denote that two functions f, g are equal up to exponentially decaying factors, in the sense that there is some constant We therefore have where S is known as the scattering matrix. The reflection and transmission coefficients r + , t + are the coefficients of the outgoing part of the field in the case u in (x) = e ik−·x , i.e. when the incident field is a plane wave from the positive x 0 direction (and reversely for r − , t − ). The following theorem was proved in [8].
Theorem 4.31. Assume that PD = D and that 0 ≤ ω ≤ K √ δ for some constant K. Let c = ∂D y 0 ψ 0 1 (t) dσ(y) and assume that c = 0. Then we have the following asymptotic expansion of the scattering matrix as δ → 0 where the error term is uniform with respect to ω.
The scattering matrix S contains two transmission peaks originating from the resonances ω α 1 and ω α 2 . Since the imaginary part of ω α 2 is very small, this corresponds to a sharp peak which will interfere with the broader peak associated to ω α 1 to create an asymmetric, Fano-type, transmission peak as illustrated in Figure 27.
θ · · · · · · u in L (a) Sketch of a screen of resonators with an incident plane wave u in . In this case, we have resonators arranged in a P-symmetric dimer that is inclined at an angle of θ to the plane of the metascreen.

Extraordinary transmission and unidirectional reflection
We now investigate further scattering phenomena of resonator screens. In particular, we will extend the analysis in the previous section (valid in the Hermitian case when v i ∈ R) to the non-Hermitian case with a PT -symmetric screen. We assume that in other words that δ 1 v 2 1 = T (δ 2 v 2 2 ). Then we can apply Theorem 3.21 to conclude that where k 0 and c are defined as in the previous section. The following theorem, proved in [11], describes the scattering behaviour of this PT -symmetrical screen of resonators.
Theorem 4.32. Assume that D is PT -symmetric, so that PD 1 = D 2 while δ 1 v 2 1 = a+ib and δ 2 v 2 2 = a−ib for a > 0, b ≥ 0. Moreover, assume that P 3 D 1 = D 2 . Let c = ∂D y 0 ψ 0 1 (t) dσ(y) and ω * = 2aC 0 11 |D1| . Assume that b|Y l | = a|c| and that ω ∈ R is in the subwavelength regime such that ω − ω * = O(δ). We then have the following asymptotic expansion of the scattering matrix: where the error term is uniform with respect to ω in a neighbourhood of ω * . In particular, we have and, at leading order, r + and r − vanish at ω = ω + and ω = ω − , respectively, which are given by Comparing Theorem 4.31 (which is valid only for b = 0) and Theorem 4.32, we see that there is an extra term in Theorem 4.32 corresponding to the non-zero gain and loss b = 0. This term is responsible for the approximate zeros of r + and r − leading to unidirectional reflection. This is numerically verified in Figure 28, where we observe that the two reflectances vanish on different sides of the critical frequency ω * . Moreover, at least formally, we see that when Im(ω α 2 ) = 0, the singularity of the scattering matrix S will not vanish in the case b = 0. In the Hermitian case, real resonances correspond to bound states in the continuum, which decouple from the far-field and therefore cannot be excited by incoming waves. In the non-Hermitian case, however, we can have real resonances with modes which are excited by incoming waves. Such resonances correspond to extraordinary transmission, where the transmitted field is greatly amplified. This amplification, which is impossible in the Hermitian case due to energy conservation, is possible due to the energy input in the non-Hermitian case [11].  Figure 28: Plot of the transmittance T± = |t±| 2 (blue) and reflectance R± = |r±| 2 (red) as functions of the frequency. Due to reciprocity we have T+ = T−. However, we observe that R+ = R−. In particular, there are points when R± vanish while R∓ is non-zero, known as unidirectional reflection. Moreover, in this non-Hermitian case the transmission coefficients are not bounded by unity and can attain large peak values, known as extraordinary transmission.

Time-modulated metamaterials
We have now seen a variety of phenomena that can occur when the material parameters vary periodically in the spatial variable x. Mathematically, we can treat the time variable t in a similar fashion. If the material parameters depend periodically on t, there can be conceptually similar phenomena which nevertheless have fundamentally different physical implications. The Helmholtz equation we have studied so far originates from the scalar wave equation when posed in the frequency domain, and is valid only when the material parameters are constant in time. To study the time-dependent case, we return to the wave equation Here, κ(x, t) and ρ(x, t) are the material parameters. We consider the case of a finite collection of resonators D = D 1 ∪ ... ∪ D N in d = 3 spatial dimensions. We assume that the modulation is only performed inside the resonators, so that We assume that κ i ∈ C 1 (R) for each i = 1, ..., N . We will consider the case when the modulation of 1/ρ and 1/κ consist of a finite Fourier series with a large number of Fourier coefficients: for some M ∈ N satisfying M = O δ −γ/2 for some 0 < γ < 1. Moreover, we assume that they are periodic in t with period T and frequency Ω = 2π T .
As before we define the (time-dependent) contrast parameter and wave speeds as for i = 1, ..., N , where δ 1. The notion of frequency is slightly altered in this time-modulated setting. Since the wave equation (4.35) is periodic in t, we can apply the Floquet transform in t and obtain the differential problem    ∂ ∂t (4.37) It is apparent that ω plays the role of the quasiperiodicity in the Floquet-Bloch theory, and we will refer to ω as a quasifrequency. Since ω is defined modulo Ω, we define the time-Brillouin zone Y * t := C/(ΩZ). Observe that we allow complex quasifrequencies; real ω correspond to solutions which are bounded in t, whereas non-real ω correspond to exponentially increasing or decaying solutions.
Due to the periodic nature of Y * t , the usual definition of subwavelength frequencies does not apply to quasifrequencies. For example, in the particular case when Ω = O(δ 1/2 ) (which will be of interest later on), the whole Brillouin zone scales as O(δ 1/2 ), meaning that all quasifrequencies tend to zero as δ → 0. In order to distinguish between these quasifrequencies, we introduce the following definition.
Definition 4.33 (Subwavelength quasifrequency). A quasifrequency ω = ω(δ) ∈ Y * t of (4.37) is said to be a subwavelength quasifrequency if there is a corresponding Bloch solution u(x, t), depending continuously on δ, which can be written as The following theorem, proved in [19], gives the capacitance matrix approximation to the subwavelength quasifrequencies as δ → 0.
Theorem 4.34. Assume that the material parameters are given by (4.36). Then, as δ → 0, the quasifrequencies ω ∈ Y * t to the wave equation (4.35) in the subwavelength regime are, to leading order, given by the quasifrequencies of the system of ordinary differential equations for y i (t), for i = 1, ..., N .
The matrix appearing in the left-hand side of (4.38) is the capacitance matrix. We can rewrite (4.38) into the following system of Hill equations  where the vector Ψ and the matrix M are defined as with W 1 , W 2 and W 3 being the diagonal matrices with diagonal entries for i = 1, ..., N . The time-modulated case shares many similarities with the non-Hermitian case (with complex parameters) studied before. Both these cases have energy input and output to the system, and not surprisingly we can find exceptional points in the time-modulated systems. Figure 29 shows the emergence of an exceptional point in a dimer structure, similar to Figure 10 but instead due to the time-modulation.

Near-zero metamaterials
In Section 4.1 we saw the emergence of exotic parameter values, namely negative effective parameters, due to the small-scale structure of the metamaterial. In the present section we will observe effective material parameters which are close to zero, in which case the Helmholtz equation reduces to the Laplace equation and wave propagation occurs without phase change (corresponding to "infinite" phase velocity). We consider a two-dimensional infinite honeycomb crystal in two dimensions depicted in Figure 30. We let the lattice Λ be generated by the lattice vectors

Near-zero refractive index in honeycomb crystals
We assume that each unit cell contains two resonators, D = D 1 ∪D 2 , such that each resonator is invariant under rotation by 2π/3 and so that D is invariant under rotation by π.
The dual lattice Λ * is generated by α 1 and α 2 given by The points in the Brillouin zone are called Dirac points. Next, we will study the band functions and Bloch modes around these points. For simplicity, we only consider the analysis around the Dirac point α * := α * 1 , the other point having a similar behaviour.  At α = α * , the generalized capacitance matrix C α has an eigenvalue of multiplicity 2: λ α * 1 = λ α * 2 . The next result shows that this asymptotic degeneracy is in fact an exact degeneracy, and moreover that the band functions intersect in a Dirac cone at this point (see Figure 31a) [14].    Next, we will see that the near-zero property follows as a direct consequence of the Dirac cone. We consider a homogenization setting close to the Dirac points. We rescale the unit cell by replacing Y with sY for some small s > 0. To fix the order of the resonant frequencies, i.e. ω α i = O(1), we assume that δ = O(s 2 ) as s → 0. We then have the following result [20].
Theorem 4.36. For frequencies ω close to the Dirac frequency ω * , namely, ω − ω * = β √ δ, the following asymptotic behaviour of the Bloch eigenfunction u α * /s+α s holds: where the macroscopic field [ũ 1 ,ũ 2 ] T := [Ae iα·x , Be iα·x ] T satisfies the two-dimensional Dirac equation The system of Dirac equations can be considered as a homogenized equation for the honeycomb structure. Eachũ j satisfy the Helmholtz equation ∆ũ j + (ω − ω * ) 2 µ 2ũ j = 0. (4.41) In particular, at ω = ω * this effective equation reduces to the Laplace equation corresponding to effective zero refractive index. Equation (4.41) describes the large-scale behaviour of the eigenmodes, illustrated in Figure 32. We emphasize that in addition to this large-scale behaviour, there will be small-scale oscillations described by the functions S α * ,k D as illustrated in Figure 31. κ j (t) = 1 + sin Ωt + 2πj 3 · · · · · · · · · · · · Figure 33: Illustration of the trimer honeycomb lattice, with phase-shifted time-modulations inside the trimers.  The Dirac cone observed in the previous section is located at the corner of the Brillouin zone, which means that corresponding modes do not lie in the radiation continuum. For physical structures of finite dimensions, these waves will be confined to the structure and will not interact with incoming waves. Physically, this can be seen as an impedance mismatch, leading to poor transmission, due to the fact that this structure corresponds to a single-near zero material, where only one of the material parameters vanish.
In order to overcome this low transmission, it is desirable to create structures with Dirac cones at the origin of the Brillouin zone. Corresponding modes lie in the radiation continuum, and will therefore interact with incoming waves. In order to achieve this, we study a time-modulated honeycomb structure as illustrated in Figure 33. Using the theory from Section 4.8 we can compute the (quasi-)band structure associated to this material. In Figure 34 we see the static (folded) band structure, with six bands in the subwavelength regime. For a particular modulation strength, we see that a Dirac cone degeneracy may appear at the origin of the Brillouin zone, enabling a double-near zero material around this point.

Concluding remarks
In this review, we have studied several Helmholtz scattering problems posed in the subwavelength regime and repeatedly encountered the concept of capacitance. We have studied the mathematical properties of the generalized capacitance matrix, in both finite and infinite, periodic settings. We were then able to use the capacitance formulation to study several different interesting physical phenomena in the field of subwavelength metamaterials. This formulation emerged from a functional analytic approach, where the generalized capacitance matrix describes the perturbation of the kernel of a non-linear integral operator for asymptotically small parameter values. Similar approaches apply to a plethora of other subwavelength scattering problems, including high-contrast dielectric particles, plasmonic nanoparticles and Helmholtz resonators [3,24,25,26,27]. Suitable capacitance formulations can thereby be used to characterize a wide range of subwavelength resonance phenomena.
Approximating classical wave systems in terms of generalized capacitance matrices shares similarities with the tight-binding approximation that is used widely in quantum theory, in the sense that both formulations provide a discrete approximation to a continuous differential problem. As we have seen, there are nevertheless fundamental differences, for instance due to the strong interactions between the subwavelength resonators. As observed in Remark 2.17, the correspondence between the capacitance formulation and the tight-binding approximation holds only in the case of dilute resonators. Even in the dilute case, long-range interactions between subwavelength resonators are relatively strong and nearest-neighbour approximations are not generally appropriate. The strength of the capacitance formulation is that the capacitance matrix accounts for these strong interactions, thereby providing a unified mathematical model to study challenging problems in subwavelength physics.

A Abstract capacitance matrix
In this section, we describe how the capacitance formulation emerges from the structure of a general integral operator A, which describes a subwavelength resonance problem. We work with the functional analytic approach described in Section 1.3. At δ = 0, we assume that the operator A(ω, 0) has a characteristic value ω = 0 of multiplicity 2N , admitting the following pole-pencil decomposition: where ker(A(0, 0)) = span{Ψ j }, ker(A * (0, 0)) = span{Φ j } and R is holomorphic for ω in a neighbourhood of 0. Moreover, we assume that A(ω, δ), for small but non-zero δ, satisfies Defining B(ω, δ) = ω 2 R(ω)L(ω, δ) yields The characteristic values are therefore determined by (A.2), which is in general a non-linear eigenvalue problem since L and B depend on ω. If we consider the subwavelength resonances, we have B = O(ω 2 δ) uniformly for ω and δ around 0. Similarly, we have L = L 0 +L, whereL = O(ωδ). Therefore, the subwavelength resonances are approximated by the eigenvalues of the finite-rank operator −KL 0 whose restriction to ker(A(0, 0)) is given by the generalized capacitance matrix: Then the characteristic values satisfy ω n = ± λ n + O(δ), where λ n are the eigenvalues of C.
As an example of this formulation, for a finite collection of N resonators in d = 3 (as considered in Section 2), the operator A is given by From this, it is straightforward to compute C ij = δiv 2 i |Di| χ ∂Di , ψ j , as defined in (2.16) in Definition 2.6.