A mathematical and numerical framework for gradient meta-surfaces built upon periodically repeating arrays of Helmholtz resonators

In this paper a mathematical model is given for the scattering of an incident wave from a surface covered with microscopic small Helmholtz resonators, which are cavities with small openings. More precisely, the surface is built upon a finite number of Helmholtz resonators in a unit cell and that unit cell is repeated periodically. To solve the scattering problem, the mathematical framework elaborated in [Ammari et al., Asympt. Anal., 2019] is used. The main result is an approximate formula for the scattered wave in terms of the lengths of the openings. Our framework provides analytic expressions for the scattering wave vector and angle and the phase-shift. It justifies the apparent absorption. Moreover, it shows that at specific lengths for the openings and a specific frequency there is an abrupt shift of the phase of the scattered wave due to the subwavelength resonances of the Helmholtz resonators. A numerically fast implementation is given to identify a region of those specific values of the openings and the frequencies.


Introduction
Surfaces covered by a microscopic structure can display unforeseen physical properties, which can be applied in the everyday life to have advantageous effects like anti-reflection coating and high-efficiency light absorbers, or enhancers. These are called meta-surfaces and have been subject to intensive research [15,24,13,31,14,22,23]. We are particularly inspired by the findings in [28,27]. In [27], the authors built a gradient meta-surface covered with microscopic small gold plates of different sizes, which are energized by an electric current, and afterwards the authors physically illuminated that surface by an incident light with a particular frequency and particular angle with respect to the ground, and emphasized a reflection with an altered outgoing angle. For certain angles of incidence, they achieved an absorption of the impinging wave. More importantly, they observed at the resonant frequencies an abrupt shift of the phase of the scattered wave. This unusual behaviour of the phase shift at muti-resonances is intrinsic to gradient meta-surfaces [11].
Here, we mainly consider acoustic waves. Our objective is to uncover the behaviour of a meta-surface built upon microscopic Helmholtz resonators, one such Helmholtz resonator can be any cavity with a small hole, as for example a water bottle with a small opening compared to its size. Such a water bottle admits a high pitched noise when blown upon, that is an effect of the resonance. To be more precise in our design, the structure is a periodically repeated ensemble of rectangular Helmholtz resonators each with a different size, where the opening is placed on the center of their ceiling. As it was the case in the paper with the gold plates, we let a plane wave hit the gradient meta-surface with a certain angle and frequency. We can expect that for appropriate frequencies and appropriate angles of incidence, the resonances of the cavities interact with each other and produce a modified reflected wave. Although we work here in two dimensions only, it is enough as pointed out in [27], because there the third dimension, which corresponds to the depth, has no measurable influence.
A solution for the scattering problem is provided in Theorem 3.1. The numerical results show an abrupt phase-shift for the scattered wave. Applying the techniques developed in [6,8,4], we provide an asymptotic formula in which the resonance is involved in the scattered wave. That equation depends on the controllable parameters, i.e., the incident frequency and the lengths of the openings of the Helmholtz resonators. It also depends on a matrix, which captures the coupling between the Helmholtz resonators. Controlling such parameters, we can exploit the resonance effect, which leads to a precise description of the unusual scattering and absorption properties of the gradient meta-surface. This is further investigated in our numerical applications at the end of this paper. From the proof of Theorem 3.1, we can extract a formula for a numerical evaluation of the scattered wave in the near-field and there we observe other remarkable outcomes like an enhancement of the amplitude of the scattered field, which is reminiscent of the results in [14]. Finally, it is worth emphasizing the connection between our results and those obtained for periodic arrays of narrow slits in the series of papers [21,20,19,18,16,17].
This paper is organized as follows. In Section 2 we prepare the mathematical foundation for the main result. There we define the geometry of the structure, the equation for the plane wave and the partial differential equation which models the solution to the problem, and make a short remark about the physical meaning for the boundary conditions for acoustic walls. In Section 3, we show the main results and comment on the appropriate choices of parameters, like the coupling matrix and the scaling factor and describe a way to utilize the resonance effect. In Section 4, we state and prove Theorem 3.1. In Section 5, we first explain how we implemented the numerics in Matlab and how we tested them. Then we present numerical applications of the results of Theorem 3.1 for some certain geometries. In Section 6, we conclude the paper with final considerations, open questions and possible future research directions.

Preliminaries
We first fix the geometry of our problem, i.e., the lengths and heights of our Helmholtz resonators, and the positions and lengths of their openings. This allows us to provide a mathematical model for the resulting wave when we imping the structure with an incident wave with a particular angle to the ground and a particular wave vector. That model is built upon a partial differential equation on a unit strip with an outgoing wave condition and a quasi-periodicity condition on both sides of the unit strips. With that in hand we then can state the main result of this paper.

Geometry of the Problem
Let δ > 0, denote a small parameter, let N ∈ N denote the number of Helmholtz resonators, then we denote a single Helmholtz resonator in the unit strip Y : We set every Helmholtz resonator D i to be a rectangle of height h i > 0 and of length l i > 0 and its center to be located at (ξ i , h i/2), where ξ i ∈ (− δ /2 + l i/2, δ /2 − l i/2). Thus the lower boundaries of the Helmholtz resonators intersect the horizontal axis, but the boundaries on the side do not intersect the side of the unit strip. Every Helmholtz resonator D i has also an opening gap Λ i := (ξ − ε i , ξ + ε i ) × {h i } located at the center of its ceiling, that is at (ξ i , h i ) and it has a radius of ε i > 0, hence a diameter of 2 ε i . This forms the geometry in our unit strip Y. See Figure 1 for an illustration of an example with N = 3.
Our gradient meta-surface in R 2 is built upon periodically repeating the unit strip along the horizontal axis. We are now interested in the result when we imping our gradient meta-surface with a harmonic incident plane wave U k 0 : R 2 → C given through where i is the imaginary unit, I 0 ∈ R denotes the intensity of the incident wave, and k 1 ∈ R and k 2 < 0 are the horizontal, respectively vertical, components of the wave vector (k 1 , k 2 ) ∈ R 2 . Let k 2 := k 2 1 + k 2 2 be the length of the wave vector. See Figure 2 for an illustration of the gradient meta-surface. The wave vector has the two components k 1 ∈ R and k 2 < 0 and the length k.

Mathematical Model for the Scattering Problem
We denote the wave function, which results from letting our geometry interacts with the incident wave, by U k : Y → C. We model U k as a solution to the following partial differential equation: where · | + denotes the limit from outside of D i and · | − denotes the limit from inside of D i , and ∂ ν denotes the outside normal derivative on ∂D i , for all i ∈ {1, . . . , N}. The first condition, that is, the Helmholtz equation, represents the time-independent wave equation, and arises from the wave equation by separation of variables in time and spatial domain. The second and third conditions represent the transmission conditions. The forth and fifth conditions represent rigid walls. It is also known as the Neumann condition. For acoustic waves, the sixth condition represents a sound cancelling ground layer. It is also known as the Dirichlet condition. Similar to diffraction problems for gratings (see, for instance, [5,9] and the references therein), the above system of equations is completed by a certain outgoing radiation condition imposed on the scattered field U k s := U k − U k 0 and a quasi-periodicity condition on U k , that is, Both conditions follow from asserting that the scattered wave U k s behaves in the same way in the far-field and in periodicity as the horizontally reflected incident wave. We remark here that the horizontally reflected incident wave, that is the function I 0 e −i k 1 x 1 e i k 2 x 2 , is the scattered solution to (2.1) in absence of Helmholtz resonators.
We further remark that in the general case, where the incident wave is a superposition of plane waves, we can decompose the incident wave using Bloch-Floquet theory, see, for instance, [25,3]. We obtain a family of problems to solve, each one with its own outgoing radiation condition. The final solution is then the superposition of all of these solutions.

Main Result
We assume that which originates from the fact that δ 2 k 2 has to be smaller than the first nonzero Laplace eigenvalue with Neumann conditions for every Helmholtz resonator D i . Furthermore, we assume that This is due to Lemma 4.1, which expresses the behaviour of certain auxiliary functions. This condition can be relaxed, but then the formulas get more complicated.
Theorem 3.1 Let the matrix R δk ∈ C N×N be defined by (4.17), let the constants r i , r ∂ i ∈ C, for i ∈ {1, . . . N}, and the constant r ex ∈ C be defined through Lemma 4.1. Let ε := max(ε 1 , . . . , ε N ) > 0 be small enough. Then for all z ∈ Y, where z 2 is large enough, we have that for a constant C > 0 independent of z, where · 2 denotes the Frobenius norm for matrices and where I s ∈ C is given through where f δk ∈ C N is a vector with the components and where |D| is a N × N diagonal matrix with diagonal entries |D i |.
The numerical implementation of the matrix R δk and the constants r i , r ∂ i , r ex is explained in Section 5.1. The value δ > 0 only appears in the form δk 1 or δk 2 in the theorem, thus we can use δ only to scale the incoming wave vector (k 1 , k 2 ). Consider that the geometry scales with δ, and it does so in both spatial dimensions simultaneously.
The value I s depends on the inversion of an N × N matrix Q(δk), and for some complex values of δk, which are close to the physical resonance values of the system, we obtain a blow-up in the entries of Q(δk) −1 . In our set-up we do not allow for a complex wave vector, only for real values for δk, thus we never encounter a blow-up of Q(δk) −1 . Nonetheless, when δk approaches the real part of those singular complex values, we see numerically a local extremum in the entries of Q(δk) −1 and hence a significant effect of the gradient meta-surface on the scattered wave. Now Q(δk) is of the form: a diagonal matrix added to R δk , and we can vary the entries in the diagonal matrix by varying all ε i . This gives us a control on the singular behaviour of Q(δk). Intuitively, to achieve a highly singular behaviour, we require all the eigenvalues of Q(δk) to be close to zero, which we obtain when the entries of the diagonal matrix 1 (δk) 2 |D| −1 + 2 π diag(log(ε i /2)) are close to the opposite of the eigenvalues of R δk . In the simplest case, R δk has one real eigenvalue λ R δk of multiplicity N, then we vary those ε i so that − 1 |D i |(δk) 2 + 2 π log(ε i /2) is equal to that eigenvalue for all i = 1, . . . , N. This leads to the equation Since R δk will generally not have only one real eigenvalue with multiplicity N, we also vary (δk) 2 in the hope to obtain an eigenvalue with multiplicity at least 2. In the numerical simulations we elaborate on how we determine those ε i on the considered geometries.
Note that the denser the eigenvalues of R δk are, the broader is the range of frequencies for which one can observe a significant effect of the gradient meta-surface on the scattered wave.

Proof of the Main Results
The outline of the proof is as follows. Before we can begin with the actual idea to prove the result, which starts in Subsection 4.4, we need a bundle of auxiliary functions, like the Neumann functions, which describe the solution to (2.1) with closed openings, and for these Neumann functions we need the fundamental solution to the Helmholtz equation, which we denote here by the Gamma-function. In Subsection 4.2, we formulate integral equations, which allow for a numeric evaluation of these Neumann functions.
Then we begin the proof by defining a scaled form of the resulting wave vector. With that microscopic view, we use Green's identity and the conditions from (2.1) to establish an integral equation on functions defined on the openings, which is solvable due to an invertible hypersingular integral operator, defined in Subsection 4.3. Its solution depends on the constant matrix R k given in Theorem 3.1. Using Green's identity again, we can recover the behaviour of the resulting wave on the far-field from its behaviour on the openings. Re-establishing the macroscopic view, we have then proven Theorem 3.1.

The Gamma-Function and the Neumann-Function
Let the wave vector satisfy the condition Then we define the quasi-periodic fundamental solution Γ k + to the Helmholtz operator + k 2 by for z, x ∈ Y, z = x. The sum in Equation (4.2) can also be defined if k does not satisfy the above condition (4.1), for this we refer to [4, Section where δ z (·) denotes the Dirac mass in z, and Γ k + is quasi-periodic in its second variable, that is, . The + in Γ k + originates from the fact that for x 2 large enough, Γ k + (z, x) is proportional to e +i k 1 x 1 for a fixed x 2 . For x close enough to z we have that Next we define the Neumann-functions N k Since D i is assumed to be a rectangle, we can express N k i (z, x) as a conditionally convergent sum for x ∈ D i . This can be readily assembled through [3, Proposition 2.7] and [12, Section 3.1]. We also see that N k i is well-defined except when k 2 is a Laplace eigenvalue with Neumann boundary conditions. We exclude these values for k from now on. We can express N k i (z, x) for x close enough to z, and k smaller than the first non-zero Laplace eigenvalue with the Neumann boundary condition through Lemma 2.9], with H 3 /2 being a Sobolev space. We know that R k i (z, ·) is analytic in k for k smaller than the first non-zero Laplace eigenvalue with Neumann conditions in the domain D i , see [3,Proposition 2.7].
The function N k where it is additionally required that the quasi-periodicity condition is satisfied, that is, N k , as well as the radiation condition, that is, We note here that the normal to ∂D i still points outside. We can write the Neumann-function N k + as where the remainder R k + exists in H 1 loc (Y), we refer to [3,Chapter 7]. An integral equation is given in the next section.
Until now we have assumed that the variable z is an element in an open domain, but we also need that z is located at the boundary. When z is at the boundary, the Dirac measure in the Helmholtz equation has its mass on the boundary as well and since an integral over the domain covers only half that mass, we have to multiply by a factor of 2 to recover the whole mass, and with that, the spatial singularity at z is twice the singularity when z is in the domain. This changes the Neumann function, as well as the remainders. We denote these variations by Note here that the singularity in k has no factor of 2. We see this by considering that the function N k i, ∂ (z, ·) − 1 π log(|z − ·|) satisfies the Helmholtz equation with a zero right-hand side with smooth Neumann boundary conditions and thus can be formulated as the convolution of N k i and the Neumann boundary data, see [3,Section 2.3.5.]. Since N k i has the singularity 1 k 2 1 |D i | so does then N k i, ∂ . This shows that the quotient between R k i and R k i, ∂ as well as the quotient between R k + and R k +, ∂ might not have a simple relation like them being equal to 2 or 1 /2 after taking a limit to the boundary.

Integral Equation for
Let us consider the integral equation for R k +, ∂ (z, x) first. Using Green's identity on ( + k 2 )N (−k 1 ,k 2 ) +,∂ (x, y) R k +, ∂ (z, y) for z, x ∈ ∪ i ∂D i , over the unit strip Y, but outside the Helmholtz resonators, we obtain with the splitting N where we used the quasi-periodicity condition, the Dirichlet ground condition and the radiation condition. Here we note that the first-named integral has R . This is due to the quasi-periodicity conditions. We can correct this, by using Green's identity once again to obtain that which then leads to the final integral equation for R k +, ∂ for z, x ∈ ∪ i ∂D i , With this at hand, we can formulate an equation for x being outside of all D i , through analogous arguments, and obtain If z is outside of all D i , the factor 2 does not appear in the splitting of the Neumann function, and we can still use the same arguments as before and obtain for R k + for z outside of all D i and x on the boundary that and if x is outside the boundary we have With these integral equations (4.6) to (4.9) and the formula for the propagating part in Γ k + , see Equation (4.2), we can readily show the following lemma: Lemma 4.1 Let k satisfy Condition (4.1). There exists a C > 0 such that for all z ∈ Λ i and for x 2 → ∞ we have that where r ∂ i does not depend on z. And for all x ∈ Λ i and z 2 → ∞ we have that where r i does not depend on z. And for all z outside of all D i and x 2 → ∞ we have that where r ex does not depend on z or x.
Here we note that if k does not satisfy Condition (4.1), then the propagating part in Γ k + has a different form, which then will alter the formulas in the above lemma as well. We do not consider that case in this paper.

The Hyper-Singular Operator L ε
Let µ represent the distributional derivative of µ. We define the following spaces and their respective norms: Let ε > 0 be small enough. Then the operator L ε : X ε → Y ε is bijective, see [26,Chapter 11], and we have that (4.14) A general closed form for (L ε ) −1 is given in [6, Proposition 3.11], from which the function (L ε ) −1 [1](t) can be determined.

The Microscopic View
We define u δk : R 2 + → C by u δk (x) := U k (δ x), where U k is the resulting wave function, and analogously u δk 0 : R 2 + → C, u δk 0 (x) := U k 0 (δ x). We readily see that u δk satisfies (2.1) with a scaled geometry and a scaled wave vector (δk 1 , δk 2 ). Furthermore, we have the following radiation condition and quasi-periodicity condition: To reduce the amount of symbols we uphold the notation for Y, D i , Λ i , l i , h i , ξ i , ε i after scaling.

Collapsing the Wave-Function to the Openings
Lemma 4.2 Let i, j ∈ {1, . . . , N}, i = j and let z ∈ Λ i . Then we have that Proof Using Green's identity with the Neumann function N δk i, ∂ inside D i , we readily obtain for z ∈ Λ i that Analogously we apply Green's identity to u δk using N δk +, ∂ on Y with the properties mentioned in Section 4.1 and obtain For an explicit demonstration we refer to [6, Proof of Proposition 3.5]. Using Equation (4.10), that is, R δk +, ∂ (z, x) = e −i δk 1 z 1 e i δk 1 x 1 e i δk 2 x 2 r ∂ i + O(e −Cx 2 ) for x 2 → ∞ and z 2 = h i , where r ∂ i does not depend on x, we can infer that lim r→+∞ Setting the first and second equations to be equal for z ∈ Λ i and applying the last equation to the second one, we prove the lemma.
We define then µ i := ∂ ν u δk | Λ i , and µ := (µ 1 , . . . , µ N ) as well as R δk : with y i = ξ i + t and z i = ξ i + τ for i ∈ {1, . . . , N}. Using the singular decomposition formulas for Γ δk + , N δk i, ∂ and N δk +, ∂ , which are written down in Section 4.1, we can rewrite the equation in Lemma 4.2 as for τ ∈ (−ε, ε) and i ∈ {1, . . . , N}, as long as δk is smaller than the first non-zero Laplace eigenvalue with Neumann boundary conditions. We define f δk i ∈ C to be which expresses the zero order term of the right-hand side with respect to ε. To simplify notation, we define f δk := ( f δk 1 , . . . , f δk N ) T , where the superscript T denotes the transpose, and we define L ε and |D| to be diagonal matrices with diagonal entries L ε i and |D i | respectively. Using Taylor series, we can decompose R δk as where R δk ∈ C N×N depends on δk and ε −ε µ := (

Solving for ∂ ν u δk on the Openings
We can rewrite Equation (4.16) in the matrix-form as (4.18) up to some error of order ε, where still µ i := ∂ ν u δk | Λ i . We can solve that equation in matrix form for ∂ ν u δk | Λ i . To this end, we consider the series expansion R δk where |D| is a square diagonal matrix with diagonal-entries |D i |. Consider that (L ε ) −1 [v](t) = v π log(ε/2) √ ε 2 −t 2 for a constant function with value v ∈ C N , where ε = (ε 1 , . . . , ε N ) and the elementary operations are defined element-wise.
After applying ε −ε on both sides of the above equation we can factor out ε −ε µ ∈ C N on the left-hand side. With the definition where log(ε/2) is a N × N diagonal matrix with entries log(ε i /2), we can solve for ε −ε µ and obtain Embedding this in Equation (4.19), we obtain Considering the error term, which we had for R δk as well as for the right-hand side for the equation in Lemma 4.2, we can readily establish an error term for ∂ ν u δk of order O N ε √ ε 2 −t 2 Q(δk) −1 2 , where we used [6, Lemma 3.14].

Remark 4.3
The values δk res ∈ C for which the operator in Equation (4.18) is not invertible are resonance values. Using the generalized Rouché theorem [3, Theorem 1.2] we see that those resonance values are of order max i (ε i ) away from those values δk app ∈ C in Equation (4.19), which yield no solution µ. We can readily give a condition to determine the δk app , that is,

Solving for u δk on the far-field
Consider that for z ∈ Y, but z not element of the closure of our Helmholtz resonators, it holds analogous to the proof of Lemma 4.2, using Green's identity, the Dirichlet condition, the quasi-periodicity condition and the radiation condition, for r → ∞. We recall that N δk + (z, x) = Γ δk + (z, x) + R δk + (z, x). Using Equation (4.11), that is, R δk + (z, x) = e −i δk 1 z 1 e i δk 2 z 2 e i δk 1 x 1 r i + O(e −C z 2 ) for z 2 → ∞ and x 2 = h i , where r i does not depend on z, and using Equation (4.12), that is, R δk + (z, x) = e −i δk 1 z 1 e i δk 1 x 1 e i δk 2 z 2 e i δk 2 x 2 r ex + O(e −C (x 2 −z 2 ) ) for x 2 → ∞, where r ex does not depend on z and x, we obtain that . Re-establishing the macroscopic view and recovering the notation for Y, D i , Λ i , l i , h i , ξ i , ε, we proved the main result.

Numerical Implementation
To carry out our numerics we rely on Matlab. To find the eigenvalues we use the in-built function eig. To find the constant matrix R δk ∈ C N×N and constant vectors r i , r ∂ i ∈ C, for i ∈ {1, . . . N}, and the constant function r ex ∈ C, we have to implement Γ k + , N k i, ∂ , N k + , N k +, ∂ using integral equations, which we will explain in the following.
Unfortunately, the series expansion for Γ k + given in Equation (4.2) is not converging fast enough and thus we have to apply a fast converging alternative representation, that is the Ewald method, which is described in [10] and which is slightly altered to fit our definition of Γ k + in [4, Section 4.1]. Here we have to mention that for k ≈ 2π n − |k 1 |, for any n ∈ N, the implementation of Γ k + does not well-behave because of an instability phenomenon, which is called the case of empty resonance. To resolve this one can use different implementations like the Barnett-Greengard method [5,Section 5.4].
The computation of N k i, ∂ , is a well-studied process, see, for instance, [3, Proposition 2.8 ] and [7]. For our implementation, we used the represen- (0) is the Hankel function of the first kind of order 0, and R k satisfies the Helmholtz equation with the boundary condition ∂ ν x R k (z, x) = −2 ∂ ν x Γ k (z, x). Using Green's identity we obtain that Applying a discretization with M nodes to the boundary, and approximating the integral with the trapezoidal rule, we obtain a linear system of equation of the form ( 1 2 I M − K ) R = f, from which we obtain a discretized version R of R k (z, ·). In this process, the problem arises that the trapezoidal rule cannot efficiently approximate singular integrals, that is with increasing number of nodes M the error in the approximation might not decay. These troublesome integrals are ε −ε ϕ(t) log(|t|)dt and T 0 ϕ(t) 1/M t 2 +(1/M) 2 dt, for a smooth enough ϕ ∈ L 2 ((−ε 1 , ε 2 )). To overcome the inaccurate approximation we use the following identities: Given now R k , we obtain R k i, ∂ by This has the disadvantage that for k → 0, R k i, ∂ should converge to a finite value, but because the implementation of R k yields a low accuracy for k → 0, we have that the above implementation of R k i, ∂ still diverges for k → 0. To approximate R k i, ∂ for k = 0, we choose some non-zero values k, for which R k i, ∂ is stable and extrapolated with a second degree polynomial in k. After some testing, we found out that min( π l i , π h i ) × [0.25, 0.75] is an interval of well-computable values for k.
The computation of N k +, ∂ , follows the same idea as for N k i, ∂ , that is we use Green's identity with Γ k + and obtain an integral equation for R k +, ∂ , which is Equation (4.6), if z and x are on the boundary. Applying a discretization with M nodes to the boundary, and approximating the integral with the trapezoidal rule, this once more leads to a linear system of equation of the from ( 1 2 I M + K ) R = f, with the discretized version R of R k +,∂ (z, ·). Fortunately, Γ k + exhibits no singular behaviour in k and moreover it has the same spatial singularity behaviour as Γ k .
If x is outside of the boundary, we apply a discretization of the integral equation (4.7) which requires the solution to the integral equation (4.6). If z is outside of the boundary we proceed analogously, by first solving a discretized version of Equation (4.8) and then applying it to Equation (4.9) if x is also outside the boundaries.
Then the computation of the matrix R δk , the vectors r i and r ∂ i and the constant r ex can be done directly by applying Lemma 4.1.
We tested our implementation using a discretized version of the Helmholtz equation and checking whether R k +, ∂ and R k + satisfy those. We could not check whether those two functions satisfy the boundary condition because the numerical implementation is unstable for values z too close to the boundary. Another way to test our implementation was to double the periodicity, that is if only one Helmholtz resonator was in the unit strip, we should get the results when we make the unit strip twice as wide. We got in all cases a very small error, which we suspect to be due to the discretization. The incident wave has the wave vector (−k cos(θ), −k sin(θ)), where k = 1.663 and θ = π/6. Our geometry involves four Helmholtz resonators. The resonator parameters are given by the following table:

Numerical Applications
We set δ to be 1, since it is only a scaling parameter.
In Figure 3 we see an illustration of the unit strip Y with an incoming wave with k = 1.663 and a wavelength of 2π k ≈ 3.778. The 300 black points on the boundary of our Helmholtz resonators are the discretization points used in the numerics.
With these values we can already determine the matrix R k ∈ C 4×4 in Theorem 3.1. We are especially interested in the dependence of its eigenvalues on k ∈ (0, π) given the angle of incidence θ. Additionally, we want to show the behaviour of the eigenvalues of R k , when θ varies. In Figure 4 we have depicted for five different angles of incidence θ the real parts of the four eigenvalues of R k as a function of k ∈ (1.5, 1.755).
For θ = π 2 (normal incidence to the meta-surface), we have one nonhorizontal line and three horizontal ones corresponding to eigenvalues with constant real parts in terms of k. For each of these three eigenvalues, the imaginary part is much smaller than the real part while the eigenvalue associated with the non-horizontal line has an imaginary part and a real part of the same order. Figure 4: Real parts of the four eigenvalues of R k for k ∈ [1.5, 1.775], corresponding to five different angles of incidence θ (in green for θ = 5 6 π, in blue for θ = 2 3 π, in black for θ = π, in red for θ = 1 3 π, and in yellow for θ = 1 6 π). Here, EV(R(k)) stands for eigenvalue of R k .
In Figure 4, we pick the parameter k ≈ 1.663 and 0.8858 ∈ Real(eig(R k )) in order to determine the ε i for i = 1, . . . , 4. As elaborated after Theorem 3.1, we take for ε 1 to ε 4 the values We note here that not all ε i are physically meaningful, for example we have that ε 1 ≈ 10 −14 , but to give an understanding of the concept, we allow ε i to hold these values. Now that we have ε 1 to ε 4 we can compute I s . Since I s is a complex number, we can decompose it into a phase θ I s ∈ (−π, π] and an amplitude |I s | > 0, as I s = |I s | exp(i θ I s ). In Figure 5 we show the result. We see that for k ≈ 1.663, that is a wavelength of approximately 3.778, and for the angle of incidence θ ≈ π we have almost no absorption, that is |I s | ≈ 1. And just below that value for k we see a full absorption, that is |I s | ≈ 0. In Figure 5 (left), we see an abrupt shift of θ I s in that neighbourhood of k.
Here we want to note that in Figure 4 we have a second accumulation of the real parts of the eigenvalues at (1.65, 0.98). For those values, we achieve a singular behaviour for Q(k) too, but the values of the vector Q(k) −1 f k are moderate, and we do not see a strong resonance effect. From the proof of Theorem 3.1, specifically from Equation (4.20), we can recover an asymptotic formula for the resulting wave U k (z) for z close to the Helmholtz resonators, which is given by where Q(δk) and f δk are as in Theorem 3.1. The field inside a Helmholtz resonator is given through Equation (4.15). For k = 1.663 and θ = π/6, we obtain the solution Real(U k (z)) shown in Figure 6, where z is located close to the Helmholtz resonator D 1 to D 4 . We also see focal spots with |Real(U k (z))| > 1 on the left and right sides of the unit strip. This outcome is similar to the findings in [14]. Figure 6: The field Real(U k (z)) is depicted, where z is close enough to the Helmholtz resonators, k = 1.663 and θ = π/6. This shows the resulting scattering when applying the same set-up as the one in Figure 3.
Finally, it is worth noticing that after simulating several generic examples, it appears that in a set-up with N Helmholtz resonators, we have N − 1 resonances, which exhibit a small imaginary part and remain at given θ constant in terms of k, and one resonance, which has a large imaginary part and whose real part at given θ increases with k.

Concluding Remarks
In this paper we have established in Theorem 3.1 that the formula U k (z) − U k 0 (z) = I s e −i k 1 z 1 e i k 2 z 2 describes the scattered field in the far-field in our geometry, up to a small error, which depends on the lengths of the openings of the Helmholtz resonators. Moreover, we have seen that I s = (v δk ) T Q(δk) −1 f δk − b, for Q(δk) ∈ C N×N and f δk ∈ C N as in Theorem 3.1, for some vector v δk ∈ C N and some scalar b ∈ C, where Q(δk) is the addition of a diagonal matrix D ε,δk and the matrix R δk . D ε,δk is given through a simple formula consisting of all ε i and δk. R δk contains the intrinsic properties of the geometry. By determining the diagonal matrix according to the eigenvalues of R δk we made the matrix Q(δk) nearly singular for some values of δk, which in turn leads to an unusual behaviour for the intensity I s of the scattered field due to the resonances of the resonators. In the numerical applications we have achieved almost full absorption and almost no absorption for wavelengths of around three units and additionally, we have obtained an abrupt phase-shift in that regime. These effects are similar to those observed experimentally for gradient meta-surfaces in [28,27].
The proof requires that the wave number δk satisfies δk < 2π − |δk 1 |. In numerical simulations, this constraint leads to δk ∈ [0, 2.5]. If we look for shorter wavelengths, we need to rewrite Γ δk + , in Equation (4.2), and Lemma 4.1. This will make the formulas more complicated, but might also yield new and unforeseen results.
Since the connection between the parameters of the geometry (ξ i , l i and h i ) and the eigenvalues of the matrix R δk is not explicit, the only way to evaluate R δk is through numerically solving integral equations. Hence we are left with trying different geometries until finding a geometry with the desired properties. One future point of interest in order to find optimal configurations might be to perform a sensitivity analysis, when we change either ξ i , l i or h i . We think that the case when two Helmholtz resonators get close to touching might reveal fascinating phenomena, as it does in the case of two nearly touching plasmonic spheres, see [29,30].
Our framework is established in two dimensions. As said in the introduction, the authors of [27] saw in their physical experiments no influence of the dimension of depth, that is the third dimension. However, a physical application of our gradient meta-surface might require a formula derived from a three dimensional Helmholtz resonator. To gather such a formula one needs a higher dimensional Green's function Γ k + , that is given in [4], and associated higher dimensional hyper-singular integral operators. An application of those is given in [8].
In our mathematical model we used the Dirichlet and Neumann boundary conditions. Materials which satisfy these conditions perfectly are not commonly available. Thus different boundary conditions might model the physical problem more accurately. For example in [2], Steklov boundary conditions are studied. Moreover, mixed boundaries as they are studied in [1] move the singularity in k of the scattered wave away from k = 0. Thus with mixed boundaries we might achieve a better control of the resonance effect. Additionally, our model has boundaries of vanishing width. To implement non-vanishing widths, one has to derive the behaviour of the scattered wave on the neck of the Helmholtz resonators.
Our results can lead to many practical applications. We are looking forward to experimental realizations of the unusual effects of gradient meta-surfaces highlighted in this paper.