Landscape of wave focusing and localization at low frequencies

High‐contrast scattering problems are special among classical wave systems as they allow for strong wave focusing and localization at low frequencies. We use an asymptotic framework to develop a landscape theory for high‐contrast systems that resonate in a subwavelength regime. Our from‐first‐principles asymptotic analysis yields a characterization in terms of the generalized capacitance matrix, giving a discrete approximation of the three‐dimensional scattering problem. We develop landscape theory for the generalized capacitance matrix and use it to predict the positions of three‐dimensional wave focusing and localization in random and non‐periodic systems of subwavelength resonators.


F I G U R E 1
We develop a landscape theory to characterize low-frequency wave scattering by systems of high-contrast subwavelength resonators.(A) A randomly distributed array of 15 spherical subwavelength resonators.(B) One of the subwavelength resonant modes of the system, whose amplitude is focused in a region of three-dimensional space by the resonators.
tendency for waves to be focused and localized by the material.In extreme cases, this can amount to the complete absence of diffusion (known as Anderson localization). 1 This behavior has been observed in a variety of systems, including electromagnetic [2][3][4] and acoustic 5,6 waves, as well as quantum systems. 7In spite of the extensive body of theory that was developed to describe this phenomenon during the second half of the 20th century, 8 a notable question remained: given a random material, is it possible to predict where wave energy will be localized?Landscape functions were developed to solve this problem. 9he theory of landscape functions provides an upper bound on the eigenmodes of a differential problem. 9When this bound is tight, the eigenmodes are forced to inherit the shape of the landscape function.Hence, if the eigenmode in question is either exponentially localized or strongly focused (in the sense that it has a significant local maximum), then the landscape function can be used to predict the possible locations where it could be localized.This idea has been used to predict localization in a range of continuous 9 and discrete 10 random systems.The link between localized modes and the landscape function has subsequently been made more precise, with 11 showing how the landscape can be used to predict the support of the localized mode.
Crucially, the constant in the upper bound given by the landscape function scales in proportion to the eigenvalue of the mode in question. 9As a result, when it comes to predicting the specific regions of localization in a random material, the method works well for the Schrödinger equation because the interesting localized modes typically occur at low frequencies.Conversely, in classical wave systems, the low-frequency modes are typically delocalized and any focusing or localization occurs at higher frequencies, meaning little useful information can be deduced from the landscape function. 12An approach for defining a high-frequency landscape was developed for discrete systems by [10].For continuous classical wave systems, an alternative solution was proposed by [12], which relies on shifting the eigenvalue of the differential operator and optimizing with respect to the energy shift, using the ideas from 11 to define localization regions.A summary of this approach can be found online at [13].
In this work, we study the scattering of scalar waves by arrays of finitely many highly contrasting material inclusions.The very large material contrast means the inclusions act as strong scatterers and exhibit subwavelength resonance.That is, resonance in a regime where the incident wavelength is much larger than the dimensions of the scatterers.As a result, these systems are special among classical wave systems as they allow for strong wave focusing and localization at low frequencies.An example of a strongly focused subwavelength resonant mode is shown in Figure 1, for a system of 15 randomly positioned spherical high-contrast material inclusions.In this case, the system does not support localization since it has open boundary conditions and energy is allowed to escape to the fair field.Nevertheless, it is clear from Figure 1 that the resonant mode's profile, when restricted to a neighborhood of the scatterers, has its amplitude strongly focused in a single region of three-dimensional space.The aim of this work is to use landscape functions to predict where this focusing will occur.As we will introduce in Section 2, the high-contrast regime considered in this work is motivated by the example of Minnaert resonance, which is the deeply subwavelength acoustic resonance of an air bubble in water. 14,15he subwavelength resonance of systems of high-contrast inclusions can be characterized using asymptotic analysis.In particular, it has been shown that the resonant frequencies and associated resonant modes of these systems are given, at leading order, by the eigenvalues and eigenvectors of the generalized capacitance matrix. 16,17This is a generalization of the capacitance matrix that often appears in electrostatics, having been introduced by Maxwell 18 to model the relationship between the distributions of potential and charge in a many-body system of conductors.The application of capacitance coefficients to the high-contrast, low-frequency setting considered here is qualitatively similar, in the sense that the generalized capacitance matrix can be thought of as capturing the strength of the coupling interactions between each of the inclusions.An alternative asymptotic strategy would be to use two-scale homogenization, which has been applied to the particular problem of localized modes in systems of high-contrast inclusions by [19, 20].However, in this work we favor the convenience of the discrete approximation provided by the generalized capacitance matrix.
We will exploit the discrete asymptotic approximation given by the generalized capacitance matrix to develop a landscape theory for predicting the positions of wave localization in these subwavelength resonant systems.We will see that, using the properties of capacitance coefficients, it is possible to prove a landscape inequality within this discrete framework.This argument is similar to the ideas of [21].We will see that this inequality still suffers from the same curse described above, whereby little useful information about an eigenmode's profile can be obtained from the landscape inequality if it has a relatively large eigenvalue.We deal with this by using a transformed version of the generalized capacitance matrix to define an "upper" landscape that captures the localization of higher frequency modes.This exploits the ideas developed by [10] and took advantage of the fact that we have a discrete approximation for the continuous system.Note, however, that although these modes occur at "higher" frequencies than the lowest frequency modes, they are still deeply subwavelength.Hence, they can still be thought of as low-frequency modes, they just lie at the upper end of the subwavelength resonant spectrum.
An important subtlety of the problem considered in this work is that it is three dimensional.Many of the previous applications of the landscape function theory have been for one-or twodimensional differential systems.One simple reason for this is that it is difficult to plot and compare three-dimensional data. 13In our high-contrast, low-frequency setting, our asymptotic model reduces the three-dimensional wave scattering problem to a discrete eigenvalue problem, with the matrix dimensions equal to the number of scatterers.Since the three-dimensional behavior is fully determined, at leading order, by the vector of values on the finite number of resonators, we can understand the resonant modes from this finite set of values.This provides a way to visualize the multi-dimensional data easily and facilitates the comparisons required to apply the landscape theory.
Another important subtlety of three-dimensional scattering problems is that it is unclear when, in general, Anderson localization will occur.For example, 22 demonstrated that there is no Anderson localization when light is scattered by a random three-dimensional ensemble of point scatterers.Comparing this with the varied results of other three-dimensional studies, such as, [23][24][25][26] highlights the subtlety and difficulty of understanding Anderson localization in three dimensions.In this work, we will study systems of finitely many resonators surrounded by free space.As discussed, the radiation conditions imposed on the far-field mean that energy is allowed to escape from the system so there is no global localization.However, we can already see from Figure 1 that this system can support strong focusing of modes.
We will begin by recalling the discrete asymptotic framework offered by the generalized capacitance matrix in Section 2. We will highlight only the details needed for this study, more comprehensive reviews can be found in [16, 17].Then, in Section 3 we will present the main theoretical results, which are landscape inequalities tailored to both the upper and lower parts of the subwavelength resonant spectrum.Finally, we will demonstrate this approach on several examples in Section 4, including both random arrangements of resonators and perturbed periodic arrays.

Problem statement
This study is concerned with the scattering of scalar waves by an array of material inclusions in three dimensions.We will denote these material inclusions as  1 ,  2 , … ,   , where  ∈ ℕ is the total number of inclusions.Since we will study the resonant properties of this system, we will refer to  1 , … ,   as (subwavelength) resonators from here on.We assume that the resonators  1 , … ,   are bounded subsets of ℝ 3 with Hölder continuous boundaries.That is, the boundary of each resonator   is such that there is some 0 <  < 1 so that   ∈  1, (meaning it is locally the graph of a differentiable function whose derivatives are Hölder continuous with exponent ).
We will use the notation  to denote the collection of all the resonators, given by the disjoint union of the  resonators We suppose that the interior of each resonator is homogeneous, such that the resonator   contains material with constant wave speed   .Similarly, we suppose that the resonators are surrounded by a homogeneous background material with wave speed .We consider the propagation of time-harmonic scalar waves with frequency  and solve a Helmholtz problem for the pressure field  ∈  1 loc (ℝ 3 ), given by where  in is the incoming wave.Here,  is the outward pointing unit normal vector,  > 0 is a nondimensional parameter, and the subscripts + and − denote the limits from outside and inside the boundary, respectively.Since  is allowed to take complex values, some care is needed to define an appropriate outgoing radiation condition.This condition is needed to specify the behavior of the solution in the far field, such that energy is only able to radiate outward and the problem is well posed.When Im  ≥ 0, the condition is the famous Sommerfeld radiation condition, which is given by In the case that Im  > 0, Equation (3) implies that () decays exponentially as || → ∞ [27, Theorem 3.6].Conversely, when Im  < 0, the resonant solutions that are of interest in this work are allowed to be exponentially growing for large ||.In this case, one natural approach is to reduce the unbounded problem to a problem posed within a large ball Ω ⊂ ℝ 3 containing all the resonators (i.e., ∪  =1   ⊂ Ω).Then, the two versions of the problem can be related using a so-called capacity operator, as outlined in [28, Section 3].With this reformulation,  can be identified using the boundary value problem on Ω.Alternatively, viewing Equation ( 2) as an eigenvalue problem of the form  =  2 , resonant frequencies are poles of the resolvent ( −  2 ) −1 .This is well defined for  2 on the real axis, thanks to the Sommerfeld radiation condition, and can be continued meromorphically into the complex plane.See, for example, 29,30 for details and examples of this process.
The parameter  in Equation ( 2) is a non-dimensional parameter that describes the contrast between the materials inside and outside the resonators.We are interested in the high-contrast limit that is characterized by  being very small.This regime is motivated by the phenomenon of Minnaert resonance, which is a strongly monopolar resonance of highly contrasting material inclusions.A canonical setting for Minnaert resonance to occur is when the material inclusions are air bubbles surrounded by water. 14,15In this case,  is the ratio of the density of air and water, having a value of approximately 10 −3 .
We will characterize the subwavelength resonance of this system (2) as an asymptotic property.First, we define a resonant mode to be a non-trivial solution  that exists when the incoming wave  in is zero.When such a solution  exists for some frequency , we say that  is the associated resonant frequency.Then, we will define a subwavelength resonant mode to be any resonant mode  which is such that its resonant frequency  depends continuously on  and satisfies This condition will allow us to perform asymptotic analysis in the high-contrast limit  → 0 in Section 2.3.We will also assume that the wave speeds are approximately constant in the sense that  = (1),   = (1), and ∕  = (1) for all  = 1, … ,  as  → 0.

Boundary integral operators
We will use boundary integral operators to represent the solutions to Equation ( 2).This has two main advantages.First, our subsequent analysis will hold for a broad collection of different shapes of resonators (provided they have sufficiently smooth boundaries).Second, by taking an appropriate ansatz in terms of boundary integral operators, we can reduce the three-dimensional differential problem (2) to a problem posed only on the boundary of the resonators .
The key boundary integral operator that we will need is the single-layer potential.This is a Green's function operator that uses the standard Helmholtz Green's function, given by Then, the single-layer potential is the operator    ∶  2 () →  1 loc (ℝ 3 ) given by Various properties of the single-layer potential and its use in scattering problems can be found in [27, 31].One key property that we will need is that  0  is invertible as a map from  2 () to  1 ().The other useful fact is that the single-layer potential can be used to represent a resonant mode , which is a solution to Equation (2) in the case that  in = 0, as where ,  ∈  2 () are density functions that need to be found.The value of this representation is that a solution of the form (7) necessarily satisfies the Helmholtz equations and the radiation condition in problem (2).Hence, it remains only to find ,  ∈  2 () such that the transmission conditions on the boundary  in Equation ( 2) are satisfied.

Asymptotic methods
Under the assumption of the asymptotic regime (4), we can make use of the well-known lowfrequency asymptotic expansions of the single-layer potential.In particular, we have that as  → 0, with convergence in the sense of the operator norm for bounded linear operators on  2 ().We will not recall all the details of the asymptotic analysis here and, instead, refer the interested reader to Section 3.2.3 of [17].The main idea is that in the limiting case, when  =  = 0, we have an operator whose spectrum is straightforward to understand.In particular, Equation (2) reduces to Laplace's equation with Neumann boundary conditions on each of the  disjoint domains  1 , … ,   .It is clear that this operator has an -dimensional eigenspace (whose eigenfunctions are constant on each of the individual resonators).Then, non-zero  > 0 and  can be studied as an asymptotic perturbation of this limiting case.In particular, we can prove the existence of  subwavelength resonant frequencies, see Theorem 3.2.4 of [17].
Lemma 1.A system of  subwavelength resonators supports exactly  subwavelength resonant frequencies with positive real parts, counted with their multiplicities (i.e., a single eigenvalue should be counted once, a double eigenvalue should be counted twice, and so on).
In the  → 0 and  → 0 limit, we have an -dimensional eigenspace spanned by the functions where Here,    is the characteristic function of   .We can use these functions to develop an asymptotic characterization of the subwavelength resonant modes.The key quantity that emerges from the asymptotic expansion (8) are capacitance coefficients, which are real numbers   given by for ,  = 1, … , .These are the standard capacitance coefficients, as introduced by Maxwell. 18o account for the resonators potentially having different sizes or containing different material parameters, we need to introduce the generalized capacitance coefficients   , given by where |  | is the volume of   .It turns out that the eigenvalues and eigenvectors of the generalized capacitance matrix determine the subwavelength resonant frequencies and subwavelength resonant modes of the system of resonators, at leading order.This is made precise by the following result, that was proved in Lemma 3.2.3 and Theorem 3.2.5 of [17].
Remark 1.A direct corollary of Lemma 2 is that if the resonators are all identical (having the same size and shape as well as the same internal material parameters), then we only need to study the capacitance matrix , as defined in Equation (11).Since  and  differ by a common  2  ∕|  | factor, their eigenvalues differ by a similar factor.As a result, the expansion for the resonant frequencies takes the form for each  = 1, … , , where  (1) , … ,  (𝑁) are the eigenvalues of the capacitance matrix (11).Crucially, however, in this special case the expansion for the resonant modes  (𝑛) in Lemma 2 is unchanged if  (1) , … ,  (𝑁) are the eigenvectors of the capacitance matrix.
We will use the asymptotic expansions presented in this section to study wave localization within a system of high-contrast resonators.In particular, the aim is to predict where the resonant modes   are localized.Thanks to Lemma 2 and Remark 1, we know that the resonant modes are determined, at leading order in the asymptotic parameter , by the eigenvectors of the (generalized) capacitance matrix.In particular, the definition of   (10) means that the functions  0  [  ]() take the value 1 on the boundary of   and vanish on the boundaries of all the other resonators.Hence, Equation ( 13) tells us that the (leading-order) values of a given resonant mode on each resonator are given by the entries of the corresponding eigenvector of the (generalized) capacitance matrix.

SUBWAVELENGTH LANDSCAPE THEORY
We wish to develop a landscape theory for predicting focusing and localization in the subwavelength regime, using the capacitance matrix approximation presented above.That is, we wish to develop upper bounds for the resonant modes, analogous to the results of [9].To do so, it is important to understand some properties of the capacitance coefficients   .Using some basic properties of the single-layer potential, which can be found in, for example, 31 we can derive an alternative representation for the capacitance coefficients, as given by the following lemma.where   is the Kronecker delta.Then, the capacitance coefficients, defined in Equation (11), are given by Given the representation from Lemma 3, it is straightforward to prove the following useful properties of the capacitance coefficients, as discussed in [32].(11), satisfy   =   and   > 0 for any ,  = 1, … , .Further, if  ≠ , then   ≤ 0 and, finally, we have that
The following property follows from the results of [33], but we include a simple proof here for completeness.Our argument is inspired by the approach of [21], who used the language of comparison matrices.

Lemma 5. The capacitance matrix is invertible and its inverse 𝐶 −1 has non-negative real entries.
Proof.The invertibility of  follows immediately from Lemma 4, since the inequality   > ∑ ≠ |  | for any  = 1, … ,  implies positive definiteness.Further, since  is positive definite we have that Now, take some  > 0 such that  > max    then we can define   as and   will have non-negative entries.As a result, we see that for any 1 ≤ , , ≤  and any  > 0. Hence, Equation ( 15) is an integral of  ×  non-negative integrands, so each entry of  −1 must be non-negative.□ With Lemma 5 in hand, we are ready to prove the main result of our landscape theory for subwavelength resonant systems.
Proof.From the fact that  =  ,1 we have for any Then, we have for any  = 1, … ,  that where the final equality follows from Lemma 5. Combining Equations ( 19) and ( 20) gives the result.□ Ultimately, when it comes to predicting the shapes of resonant modes, the inequality in Theorem 1 suffers from the same problem as described in the introduction: that when  is large the inequality (18) does not reveal any useful information about the shape of the eigenvector .This means that resonant modes belonging to the lowest part of the spectrum can be captured effectively but modes at the top of the subwavelength spectrum are likely to be described in less detail.However, we can take inspiration from the matrix   that was defined in Equation ( 16) to develop a specific landscape theory for the modes at the top of the subwavelength spectrum.To do so, we will need to introduce the notion of a comparison matrix, as used previously by, for example, 21,33 For an  ×  matrix , the associated comparison matrix is the  ×  matrix  ‡ given by Using the comparison matrix of   , for  sufficiently large, we will be able to define an "upper" landscape vector û that provides a tighter bound on the eigenvalues at the top of the subwavelength spectrum.This is a variant of the high-energy landscape developed by [10] for tridiagonal discrete systems.
Proof.We will use the notation   =  − , as was used previously in Equation ( 16).The choice of  in Equation ( 22) means we can use Lemma 4 to see that Hence, we know that   is positive definite and invertible.We now need to compare   with its associated comparison matrix  ‡  .The general result here is known as the comparison theorem, which says that for an  ×  matrix , if its comparison matrix  ‡ is positive definite, then it holds that |( −1 )  | ≤ ( ‡ )  , for all 1 ≤ ,  ≤ .This was first proved in [33] (in German).An alternative proof (in English) can be found in the Supporting information of [21].In our case,  ‡  inherits the property of being diagonally dominated from Equation ( 24), so must be positive definite.Thus, we have that Note also that if an eigenvector satisfies  = , then we have that    = ( − ) = ( − ), where the choice of  as Equation ( 22) implies that  −  > 0. As a result, we can show that for any  = 1, … , . □ The value of the upper landscape developed in Theorem 2 is that the constant ( − ) in the inequality ( 23) is now such that it is smallest when the eigenvalue  is largest, meaning that this will give the tightest bound on the modes at the top of the subwavelength spectrum.We will demonstrate this with some numerical examples in the following section.This result extends the work of [10] and is qualitatively similar to the approach taken in [12], where they similarly used an eigenvalue shift to develop an "optimized" landscape function (where the optimization was performed over the energy shift parameter).Their approach relied on using Webster's transformation 34 to transform the Helmholtz equation to a Schrödinger problem, for which energy levels can be shifted freely.

EXAMPLES
One of the practical challenges that makes it difficult to study wave focusing and localization in three dimensions is the simple fact that it is hard to visualize three-dimensional functions.This is another challenge that is overcome through the use of a discrete asymptotic approximation in this work (based on boundary integral operators and the capacitance matrix), since the three-dimensional resonant modes () are fully determined (at leading order in the asymptotic parameter ) by the values of the potentials on the  resonators.That is, it is sufficient to study the eigenvectors  (1) , … ,  () ∈ ℝ  of the capacitance matrix.These vectors can be easily plotted against their index.An example of this is shown in Figure 2. Figure 2A shows a cross-section of a resonant mode that appears to be strongly focused and has been plotted in a plane on which The focusing of a resonant mode can be characterized by the values of the field on each resonator.A system of 15 identical spherical resonators with radius 1 is modeled using the capacitance matrix.
The positions are chosen to be normally distributed, with the constraints that the resonators must not overlap and their centres must lie on the  = 0 plane.(A) The 15 th resonant mode  (15) (, , ) in the  = 0 plane.(B) The 15 th mode plotted on each of the 15 resonators (with the value given by the eigenvector of the capacitance matrix).
(C) The entries of the 15 th eigenvector  (15) of the capacitance matrix plotted on a non-physical axis.The color scale is the same for all three plots.
the centres of the resonators all lie (by design, see Section 4.1).Figure 2B shows the same mode plotted only on the resonators themselves.These leading-order constant values are given by the entries of the corresponding eigenvector of the capacitance matrix.Finally, Figure 2C shows the same eigenvector with the entries plotted against their index.Although the horizontal axis is nonphysical is this case, the plot still contains meaningful information about the extent to which the resonant mode is focused in a region of space (as opposed to its amplitude being distributed across the resonator system).The numerical computations are performed using multipole expansions.We consider spherical resonators, in which case it is convenient to express the density functions   from Equation (10)  in terms of a basis consisting of spherical harmonics.Calculating the expression of   within this basis relies on understanding the image of each basis function under the operator  0  .This is explained in detail in the Appendices of [35].Given an expression for   in terms of this basis, it is straightforward to integrate over the boundaries of the resonators to calculate the capacitance coefficients.
It is important that we consider examples where the resonators are not too close together (and certainly do not touch), as it is known that the field blows up between close-to-touching subwavelength resonators. 36,37This is not a problem for the asymptotic approximation in terms of the generalized capacitance matrix (in fact, the behavior of the capacitance coefficients of close-to-touching spheres is well known 38 ), nor does it present issues for the landscape theory.In fact, it is likely to be a mechanism for strongly focused fields as it allows the otherwise monopolar scatterers to exhibit dipolar behavior (cf. the mechanism demonstrated in [39] for elastic plates).However, close-to-touching resonators pose challenges for numerical computations, as unfeasibly large numbers of basis functions would be needed to handle the problem.

Random systems
A first example is shown in Figure 2, in which 15 identical spherical resonators are distributed randomly on a plane.Since the resonators are identical, it is sufficient to study the eigenvectors of the The shapes of the eigenvectors of the capacitance matrix can be predicted by the landscape and the upper landscape.The 15 eigenvectors  (1) , … ,  (𝑁) corresponding to a random arrangement of 15 subwavelength resonators are shown, in the same arrangement as was depicted in Figure 2. In each plot, the absolute value of the eigenvector is shown, normalized to have maximum equal to 1.The upper bounds given by the landscape  (𝑛)  and the upper landscape ( −  (𝑛) ) û are also shown in each case.
capacitance matrix (see Remark 1).The resonators all have unit radius and their centers that are chosen to lie on the plane  = 0.The  and  coordinates were drawn from a normal distribution with mean zero and variance 100.The distribution of the resonators is shown in Figure 2B.In this system, the 15 th subwavelength resonant mode is strongly focused.Figure 2A shows this resonant mode in the  = 0 plane.The color of the circles in Figure 2B indicates the value of the mode on each resonator.Finally, Figure 2C shows the corresponding eigenvector of the capacitance matrix.
Although the numbering of the resonators is completely arbitrary, the focusing in a neighborhood of resonators 3, 7, and 10 is clear (and would similarly be so if a different numbering scheme was chosen).The eigenvectors of the same system are shown in Figure 3, where all 15 eigenvectors of the capacitance matrix are compared to the landsapce and the upper landscape.We can see that the landscape (from Theorem 1) accurately reproduces the profile of the first few eigenvectors but struggles to give any meaningful information beyond this.Conversely, the modes at the top of the subwavelength spectrum can be seen to be strongly focused and their shape is predicted well by the upper landscape (from Theorem 2).
We can also use this approach to study examples, where the resonators are free to take random positions in three-dimensional space.For example, let us consider a system of 20 identical spherical resonators (with unit radius) that are distributed uniformly within a 20 × 20 × 20 cube.In this case, since the resonators are not constrained to a plane, as in the previous example, we elect not to plot the full resonant mode  (𝑛) and directly opt to plot the eigenvectors of the capacitance matrix on a non-physical axis.Five of the eigenvectors for a realization of this system are shown in Figure 4. We choose the first two modes, for which the landscape from Theorem 1 predicts the shape, as well as the top three modes, for which the upper landscape from Theorem 2 is more revealing.

F I G U R E 4
The discrete asymptotic approximation used in this work means the focusing and localization of resonant modes in three-dimensional systems can be characterized concisely.Here, 20 identical spherical resonators are positioned with centers drawn from a uniform distribution on a 20 × 20 × 20 cube.Five of the eigenvectors  (𝑛) are chosen to be plotted, with each being shown alongside both the landscape  (𝑛)  and the upper landscape ( −  () ) û.
Remark 2. It is well known that the localization of eigenvectors of banded random matrices is related to their bandwidth.In particular, eigenvectors are typically localized when the bandwidth is less than  1∕2− and delocalized when it is greater than  1∕2+ . 40With this in mind, it is interesting to note that while the capacitance matrix is not banded, its entries decay approximately in proportion to | − | −1 .In particular, [35, Lemma 4.3] shows that   is approximately inversely proportional to the distance between the centes of resonators  and .While this means that the capacitance matrix is somewhat reminiscent of a banded matrix, this decay is not fast enough for it to behave like a banded matrix with small bandwidth.This is one explanation for why some of the eigenvectors shown in Figures 3 and 4 are delocalized.Furthermore, it has been observed that truncating the capacitance matrix to give a banded matrix alters the fundamental physics of the model and gives a poor approximation of the physical system (see, e.g., [35, Section 4.2.2.]).

Periodic systems with defects
Typically, periodic systems (including periodic arrangements of subwavelength resonators) do not give rise to any interesting wave focusing or localization effects. 41Instead, they typically support purely propagating modes (often known as Bloch modes) with any insulating properties arising due to spectral gaps.However, if perturbations are made to a periodic material, then this will often lead to the creation of resonant modes that are localized to region of the defect or interface that has been created.This is a popular strategy for designing waveguides and related problems have been studied extensively.This includes systems of subwavelength resonators; several different examples of this are surveyed in [16], and references therein.
A simple example of a way to perturb a periodic system to create a localized resonant mode is to perturb the material parameters on of the resonators.This will cause one the subwavelength resonant modes to be localized in a region of the perturbed resonator.In this case, we will need to use the generalized capacitance matrix to capture the behavior of the system.This setting was studied in detail in [42], where a formula was derived for the eigenvalue of the generalized capacitance matrix corresponding to the localized eigenvector.It shows that the localized is perturbed away from the spectral band of propagating Bloch modes.This localized mode will exist for any positive perturbation of the material parameter, with the exponential localization becoming increasingly strong as the perturbation increases.
An example of this system is shown in Figure 5, where we have studied identically sized spherical resonators (with unit radius) distributed on a square lattice.Within the central resonator, the wave speed   is chosen to be double than in the other resonators.This perturbation is chosen such to be sufficiently large that the mode decays observably within the length of the finite sized array depicted here.In this case, the localized mode appears at the top of the subwavelength spectrum.Hence, we a priori expect the upper landscape to be most useful for predicting its shape.The eigenvector is shown as it is distributed in space over the resonators in Figure 5A and it is shown alongside both the landscape and the upper landscape in Figure 5B.Clearly, the landscape is not useful, so a zoomed plot showing just the mode and the upper landscape can be found in Figure 5C.The upper landscape unambiguously predicts the localization on the central resonator.

CONCLUDING REMARKS
We have shown that the landscape theory can be applied directly to the generalized capacitance matrix, which is used here as a discrete asymptotic model for the subwavelength resonance of high-contrast heterogeneous systems.Using the properties of the capacitance matrix, we were able to obtain simple direct proofs of the landscape inequality (Theorem 1).We also used a transformed version of the capacitance matrix to obtain an "upper" landscape, that can describe the top of the subwavelength spectrum (Theorem 2).We demonstrated this approach on several examples, including both random arrangements of resonators and perturbed periodic arrays.The discrete asymptotic approximation at the center of this work means that three-dimensional focusing and localization problems could be handled and visualized concisely.The value of the landscape theory is particularly apparent when dealing with large systems of many resonators.In this case, substantial computational expense would be required to compute all the eigenvectors of the capacitance matrix and check where they are localized or focused (if at all).Conversely, computing the landscape and upper landscape requires just two matrix equations to be solved and gives direct insight into where energy is likely to be focused within the system.
The two versions of the landscape inequality presented in this work mean that the profile of resonant modes occurring at either end of the subwavelength spectrum can be captured.However, this theory is still less revealing for resonant modes whose resonant frequency falls in the center of the subwavelength spectrum, especially for very large systems of many resonators.There are several notable examples for which the important localized modes appear in the middle of the subwavelength spectrum, such as the topologically protected modes that are widely used in waveguide design. 35Predicting the localization of waves in quasicrystalline wave systems is a similarly challenging problem.While landscape functions have been applied to quasicrystalline Hamiltonian systems, 43 the localized modes typically appear in the center of the spectrum (within one of the fractal collection of spectral gaps 44 ), so suffer from the issue of landscape bounds not being tight enough to reveal useful information.
While understanding wave propagation in random and perturbed periodic media is interesting in its own right and has significant applications (such as to imaging, see, e.g., 45 as well as waveguide design), it also has deep implications for the field of metamaterials and wave guides.When any micro-structured wave control device is manufactured, it will inevitably have some random errors and imperfections introduced.While small random perturbations have minimal effect on subwavelength devices, 46 when these imperfections are sufficiently large it is possible that they will hinder the device's function.It is not yet clear how landscape theory can be used to understand this breakdown in function and, by extension, be used to design devices that are robust with respect to imperfections.However, this is a promising direction for future work.

C O N F L I C T O F I N T E R E S T S TAT E M E N T
The authors declare no conflicts of interest.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The code used to produce the numerical examples presented in this work is available for download at https://doi.org/10.5281/zenodo.8134312.No other datasets were generated during this study.

F
I G U R E 5 A localized resonant mode can be created by introducing a defect to a periodic arrangement of resonators.25 spherical resonators are arranged in a square lattice and the wave speed on the central resonator is perturbed to be double that on the other resonators.(A) The value of the mode on each of the 25 resonators.(B) The profile of the mode compared to both the landscape  ()  and the upper landscape ( −  () ) û. (C) The mode compared to just the upper landscape ( −  () ) û, which gives a more insightful prediction of the mode's profile (as would be expected, for this high-frequency subwavelength resonant mode).