Absorption characteristics of large acoustic metasurfaces

Metasurfaces formed of arrays of subwavelength resonators are often tuned to ‘critically couple’ with incident radiation, so that at resonance dissipative and radiative damping are balanced and absorption is maximized. Such design criteria are typically derived assuming an infinite metasurface, whereas the absorption characteristics of finite metasurfaces, even very large ones, can be markedly different in certain frequency intervals. This is due to the excitation of surface waves, intrinsic to resonant metasurfaces and especially meta-resonances, namely collective resonances where the surface waves form standing-wave patterns over the planar metasurface domain. We illustrate this issue using a detailed model of a Helmholtz-type acoustic metasurface formed of cavity-neck pairs embedded into a rigid substrate, with geometric and dissipation effects included from first principles (R. Brandão and O. Schnitzer, Wave Motion, 97, 102583, 2020). This article is part of the theme issue ‘Wave generation and transmission in multi-scale complex media and structured metamaterials (part 1)’.

Metasurfaces formed of arrays of subwavelength resonators are often tuned to 'critically couple' with incident radiation, so that at resonance dissipative and radiative damping are balanced and absorption is maximized. Such design criteria are typically derived assuming an infinite metasurface, whereas the absorption characteristics of finite metasurfaces, even very large ones, can be markedly different in certain frequency intervals. This is due to the excitation of surface waves, intrinsic to resonant metasurfaces and especially meta-resonances, namely collective resonances where the surface waves form standingwave patterns over the planar metasurface domain. We illustrate this issue using a detailed model of a Helmholtz-type acoustic metasurface formed of cavity-neck pairs embedded into a rigid substrate, with geometric and dissipation effects included from first principles (R. Brandão  large, metasurface absorbers. To this end, we shall use a reduced multiple-scattering model of a finite-sized acoustic Helmholtz-type metasurface formed of cavity-neck pairs which are embedded into a flat rigid substrate. We have recently derived this model using matched asymptotic expansions [22], first considering the acoustic response of a single cavity-neck resonator and then generalizing to an arbitrary planar distribution of resonators based on a Foldytype multiple-scattering approximation [23,24]. The geometric details of the three-dimensional cavities and necks, which can be quite general, are encoded in lumped parameters, which are defined systematically in terms of boundary-value problems and for which we have obtained numerical values and analytical approximations [22,25]. The model also allows for an arbitrary geometric arrangement of non-identical resonators, though here we will focus on a square lattice of identical resonators. Most importantly, and in contrast to preceding models of Helmholtz resonators, the model derived in [22] includes dissipative effects from first principles; these are shown to be dominated by the viscous boundary layers in the necks of the resonators. While Helmholtz-type metasurfaces are relatively basic compared to more sophisticated concepts [26,27], they remain very common in applications: essentially similar surfaces are used as absorbing panels for indoor acoustics and as acoustic liners for jet engines; the soda-can metasurface [18] already mentioned provides a particularly simple realization. The paper is structured as follows. In §2, we formulate a multiple-scattering model of the Helmholtz-type acoustic metasurface based on the theory developed in [22,25], as well as derive a homogenization approximation of that model which is valid in some cases. In §3, we derive explicit critical-coupling conditions in the two extreme cases of a single resonator and an infinite metasurface. In §4, we study the surface waves supported by an infinite metasurface and their attenuation. In §5, we numerically explore the absorption characteristics of large metasurfaces, comparing both 'tuned' (critically coupled) and 'detuned' (low-loss) metasurfaces to their infinite counterparts and attempt to qualitatively explain some of the observed characteristics. In §6, we give concluding remarks and propose directions for further study.

Helmholtz-type metasurface (a) Geometry and physical assumptions
We consider a model acoustic metasurface formed of N identical cavities which are embedded into a flat substrate and arranged in a square lattice. Each cavity (volume l 3 ) is connected to the half-space exterior to the substrate by a cylindrical neck (radius l, length 2 hl) whose axis is perpendicular to the substrate plane. The shape of the cavity is arbitrary, except that its boundary is assumed flat in the close vicinity of the neck, which does not protrude into the cavity or exterior half-space. The spacing between the necks is denoted A 1/2 l, such that Al 2 is the area of a unit cell of the lattice. The fluid filling the cavities, neck and the exterior half-space is assumed to be a viscous and thermally conducting gas (speed of sound c, density ρ, kinematic viscosity ν). The substrate is assumed rigid, isothermal, with the fluid velocity satisfying a no-slip boundary condition. Accordingly, there is an additional characteristic length scale, l v = √ ν/ω H , corresponding to the width of the viscous boundary layer, ω H being a characteristic angular frequency to be specified below; the characteristic width of the thermal boundary layer is assumed comparable, as is the case for air. Associated with l v is the ratio δ = l v /l, which serves as a measure of material loss. We adopt a dimensionless convention where lengths are normalized by l. The geometry is shown in figure 1 for the case of cubic cavities. Relative to an arbitrary fixed origin, we denote the dimensionless position vector by r and position vectors of the neck centres in the substrate plane by r n , where n = 1, 2, . . . , N. Furthermore, we denote the vertical coordinate measured from the substrate plane by z. The pressure field will be sought as the real part of p ∞ p(r) exp(−iωt), where p ∞ is a reference magnitude, ω is the angular frequency, t is time and the reduced pressure field p(r) is a dimensionless and complex-valued phasor field.
Our interest is in studying the sound-absorption characteristics of the metasurface defined above for 1

(b) Multiple-scattering model
In [22], we first studied the limit 1 in the case of a single cavity-neck resonator, with Ω 1 and δ , the symbol henceforth standing for 'of asymptotic order'. The condition δ , namely that the viscous and thermal boundary layers are thin relative to the neck radius, was shown to be equivalent to saying that dissipation is sufficiently weak such that the cavityneck pair exhibits asymptotically significant resonance. Furthermore, dissipation is then only important in a vicinity of the resonance frequency, with viscous effects dominant over thermal effects and contributed mainly by the viscous boundary layers close to the neck. By analysing a series of distinguished limits in the δ-Ω parameter space, representing different levels of loss and frequencies increasingly close to resonance, we systematically derived a 'unified' asymptotic model for the acoustic response of a single cavity-neck resonator that holds to leading asymptotic order throughout the regime Ω 1 and δ . We then presented an intuitive generalization of this unified model to an arbitrary number and arrangement of resonators, based on a Foldy-type multiple-scattering approximation [23,24] and assuming A 1/2 , namely that the separation between the necks is large relative to their radius. This includes the case A 1, where the separation is subwavelength and comparable in size to the cavity.
In what follows, we formulate a multiple-scattering model based on the theory developed in [22,25]. As → 0, the neck openings in the substrate plane shrink to the points r n and the reduced pressure p satisfies the reduced-wave equation where f (r) accounts for a possible sound sources in the exterior half-space. On the substrate plane, p satisfies the Neumann boundary condition ∂p ∂z = 0 at z = 0 (r = r n for n = 1, 2, . . . , N), (2.2) whereas, as r → r n in the exterior half-space, p approaches monopolar singularities whose amplitudes satisfy point constraints specified below. The scattered field, p(r) − p (i) (r), satisfies an outward radiation condition as z → ∞, wherein the incident field p (i) (r) is the reduced pressure in the absence of the substrate, associated with incoming far-field radiation and the source distribution f (r). In light of the above, the reduced pressure is approximated as where α n are diffraction coefficients associated with spherical waves emitted from the necks and we introduce the reflected field p (r) (r) defined as the solution for the reduced pressure in the absence of the cavity-neck resonators. The reflected field can readily be found for an arbitrary incident field using the method of images. We shall refer to the sum p (a) (r) = p (i) (r) + p (r) (r) as the ambient field. In our examples, we consider two cases: (i) normally incident plane wave, in which case p (i) (r) = e −i 1/2 Ωz and p (a) (r) = 2 cos( 1/2 Ωz); and (ii) a monopole source located at r = r s , in which case p (i) is the free-space Green's function G(r, r s ) = |r − r s | −1 e i 1/2 Ω|r−r s | and p (a) (r) = G(r; r s ) + G(r; r s ), r s being the reflection of r s with respect to the substrate plane.
The nth diffraction coefficient is determined by requiring that the resonator associated with it reacts as if it was isolated and exposed to the effective ambient pressure field i.e. the sum of the ambient pressure and the spherical waves diffracted from all of the other resonators. As already mentioned, the response of a single embedded resonator was asymptotically analysed in [22]. Thus, consistently with the standard picture of a Helmholtz resonator, the pressure in each cavity is approximately uniform. Furthermore, the cavity pressures, denoted p n , are proportional to their associated diffraction coefficients according to α n = Ω 2 2π p n , n = 1, 2, . . . , N. (2.5) This relation physically represents compression and expansion of the gas in the cavities. In the single-resonator theory, the cavity pressure is found to be proportional to the ambient pressure. Thus, for N resonators, the cavity pressures are calculated by solving the system of equations p m e i 1/2 Ω|r m −r n | |r m − r n | , n = 1, 2, . . . , N, (2.6) where the factor of proportionality comes from the asymptotic single-resonator theory. Here,Ω is the resonance frequencȳ respectively, K, Θ and σ are lumped parameters, and¯ r and¯ v are, respectively, r and v evaluated at Ω =Ω.
The parameters K and Θ are functions of the aspect ratio h of the neck. In particular, K is the so-called Rayleigh conductivity of the neck normalized by its radius [28], which is determined by solving a potential-flow problem involving the neck geometry. The parameter Θ, which captures the role of the neck geometry on viscous resistance, is calculated as a quadrature of that same potential flow, representing the displacement of that flow by the viscous boundary layer. In [25], we derived accurate asymptotic approximations for K and Θ in the limits of small and large h, corroborating and significantly extending previous approximations in the literature; we also provided numerical results for arbitrary h. The parameter σ is a cavity shape factor whose calculation entails solving a Poisson-type boundary value problem defined over the cavity domain with the neck opening reduced to a point; in [22], this problem is defined and solved for several families of cavity shapes. Expressions and numerical results for K, Θ and σ are provided in the electronic supplementary material. 1 The quantity of interest for us is the power dissipated by the metasurface, or the averaged dissipation per resonator which can be thought of as an absorption efficiency. We shall normalize power by p 2 ∞ l 2 A/2ρc, the incident power per unit area of the metasurface in the case of a normally incident plane wave. Once the multiple-scattering problem is solved for a given incident field, the total dimensionless dissipation, denoted D, is provided as Formula (2.10) is derived in electronic supplementary material (see footnote 1) in two independent ways. First, it is shown that within the framework of the multiple-scattering model (2.10) gives exactly the rate of energy leakage through the limiting neck positions r n . Second, it is shown that (2.10) asymptotically agrees, to leading order both near and away from resonance, with a direct integration of the viscous dissipation in the boundary layers near the necks. The latter calculation is carried out by relating the cavity pressures p n in the multiple-scattering model to the matched asymptotic expansions developed in [22,25] to describe the cavity, neck and exterior regions.
We note that the multiple-scattering model herein slightly differs from that in [22]. The difference is that in [22] the frequency dependence of r and v is ignored, with Ω in (2.9a,b) approximated by the leading-order resonance frequency K 1/2 . This results in a mathematically more explicit formulation and the difference is asymptotically negligible for the purpose of calculating pressure to leading order both near and away from resonance, as well as dissipation near resonance. Nonetheless, we choose to introduce this frequency dependence here as the analysis in electronic supplementary material (see footnote 1) shows that this extends the asymptotic validity of (2.10) to off-resonance frequencies and ensures that the multiple-scattering theory and (2.10) are exactly self-consistent both physically and numerically.
For the numerical examples in this paper, we shall assume a subwavelength unit-cell area A = 1, neck parameters = 0.02 and h = 2, for which K ≈ 0.56 and Θ ≈ 3.71, and a cubic cavity for which σ ≈ 0.2874. We will vary the number of resonators N and use δ to tune the material loss.

(c) Homogenized model
The multiple-scattering model can in some cases be approximated by a homogenized model where the metasurface is represented by an effective impedance condition. A necessary, but as we shall see not sufficient, condition is subwavelength spacing, i.e. A 1/2 1/ 1/2 . Let P be a 'macroscale' exterior pressure, assumed to vary on some long length scale L A 1/2 , which may or may not be comparable to the order 1/ 1/2 free-space wavelength. On this length scale, the substrate plane z = 0 decomposes into a metasurface patch M (a square of area A = NA) and the remaining domain, say S. The homogenized problem consists of the reduced-wave equation and an outward radiation condition on the scattered field P − P (i) , where P (i) = p (i) from the multiple-scattering model. The averaged dissipation per resonator is calculated as where the integrand is evaluated at z = 0 and dA is an infinitesimal area element. In the electronic supplementary material (see footnote 1), we describe an elementary boundary element method (BEM) that we use to solve the homogenized model in the case of a finite metasurface. The homogenized model can be formally derived starting from the multiple-scattering formulation by applying the method of multiple scales in an appropriate limit process. Here, we suffice with a heuristic justification. First, we assume that the effective ambient field p (a) n experienced by the nth resonator can be approximated by the macroscale field P, so that p n = gP(r n ).
( 2.15) This implies that p n and α n are slowly varying over the lattice of neck positions. The impedance condition (2.12) then follows by comparing the normal macroscale 'flux' A ∂P/∂z at z = 0 over an area A with the corresponding microscale flux −2πα n (A /A), and using (2.5). The linear dimension of the area A is taken to be large relative to the spacing A 1/2 yet small relative to the long scale L. The validity of this flux balance hinges upon the assumed slow variation of P(r), and hence p n and α n , as well as the assumption that the spacing is subwavelength. Expression (2.14) for the averaged dissipation readily follows from (2.10) upon using (2.15). We shall see that the long-scale assumption can fail at near-resonance frequencies for sufficiently low loss, owing to the excitation of short-wavelength surface waves, whereby the homogenized model loses validity.

Critical coupling (a) Single resonator
In the case of a single resonator, the cavity pressure is p 1 = gp (a) (r 1 ). From (2.10), we find using (2.7)-(2.9a,b) that the dissipation at resonance is For a fixed geometry, dissipation at resonance is optimized by tuning the loss parameter δ such that¯ v =¯ r . We shall refer to this condition, which corresponds to a balance between dissipation and radiation damping as critical coupling for a single resonator. The corresponding explicit relation obtained in [22] is recovered by consistently approximatingΩ ∼ K 1/2 . From (2.9a,b), critically coupling a single resonator demands very low loss, δ 5/2 , resulting in a narrow resonance interval, Ω 3/2 .

(b) Infinite metasurface
Consider next the scattering problem where an infinite metasurface is exposed to a normally incident plane wave. Strictly speaking, scattering problems involving infinite metasurfaces can be ill-posed in frequency intervals where surface waves are supported, that is without generalizing the radiation condition to account for this. For a dissipative metasurface, it seems sufficient to demand that the solution remains bounded. Under this assumption, symmetry of the metasurface and incident field implies that the cavity pressures in the multiple-scattering model are all equal, while the macroscale pressure in the homogenized model is a function of z alone. Let us first assume the homogenized model. In light of the above, the macroscale pressure can be written where the first term is the incident plane wave and the second term is a reflected plane wave, R being a reflection coefficient. The reflected plane wave should not be confused with the reflected field p (r) , which in this case is the reflected wave with R = 1. The deviation of R from unity represents the aggregated effect of the spherical waves emitted from the necks. Indeed, the effective impedance condition (2.12) gives the reflection coefficient R as Defining the collective radiation-damping factor and noting that r R , the reflection coefficient (3.3) can be approximated as In [22], we derive an expression for R which is asymptotically equivalent to (3.5), assuming only that A 1 as → 0, i.e. that the spacing is subwavelength, specifically comparable to the characteristic cavity size. (Minor typos in [22] are corrected in electronic supplementary material; see footnote 1). Rather than relying on a homogenized model, the equivalent result in [22] is derived directly from the multiple-scattering model. The derivation entails the asymptotic approximation of a lattice sum. Comparison of these two distinct derivations confirms that, in the present special scenario of a plane wave normally incident on an infinite metasurface, the homogenization approximation is valid for subwavelength metasurfaces regardless of the level of loss. We will see that this is not the case if either the metasurface is finite or the incident field is not a plane wave.
According to (3.5), R can be made to vanish at resonance by tuning the loss parameter δ such that¯ v =¯ R , where¯ R is R evaluated at Ω =Ω. This is the condition for critical coupling of an infinite metasurface exposed to a normally incident plane wave. The corresponding explicit condition obtained in [22] is recovered by consistently approximatingΩ ∼ K 1/2 . For a different perspective, consider the dissipation per resonator of the infinite metasurface, denoted D ∞ , which can be calculated in the homogenization approximation using (2.14) and (3.3). At resonance, D ∞ is well approximated byD which can be optimized for a fixed geometry by tuning δ to satisfy¯ v =¯ R at resonance, givinḡ D ∞ = 1 as expected for perfect absorption and given our normalization of the dissipation. As in the single-resonator case, the present critical coupling condition also corresponds to a balance between radiative and dissipative damping. The radiative damping in the metasurface case, however, is enhanced by an order 1/ factor relative to the single-resonator case. Accordingly, critically coupling an infinite metasurface requires more loss, δ 3/2 , resulting in a wider resonance, Ω 1/2 . We also see from (3.5)

Surface waves
We continue to consider the case of an infinite metasurface, now looking for homogeneous surface-wave solutions that exist for real frequencies in the absence of an incident field and decay exponentially away from the substrate. Such solutions satisfy quasi-periodicity, i.e. that p exp{−iκ · r} is periodic over the two-dimensional lattice whose vertices are the limiting neck positions r n , wherein κ = κ xêx + κ yêy is a Bloch wavevector parallel to the substrate plane with components κ x and κ y in the directionsê x andê y parallel to the lattice vectors. For a prescribed real frequency Ω, the Bloch wavevector κ is an eigenvector which for δ > 0 is expected to be complex valued, meaning that the surface waves exponentially attenuate parallel to the substrate plane. While defined for an infinite metasurface, we shall refer to these homogeneous solutions in the next section when interpreting the absorption characteristics of large metasurfaces. We first assume the homogenized model. In that case, the dispersion relation is isotropic, meaning that we can write κ = κî, κ being a complex Bloch wavenumber andî an arbitrary direction in the substrate plane. The reduced-wave equation (2.11), quasi-periodicity and the condition that the solution decays away from the substrate together imply that the macroscale pressure possesses the form with the decay condition requiring −1 b = Re(κ 2 − Ω 2 ) 1/2 > 0, b being the length scale of decay normal to the substrate. We also define the propagation length p = 1/Im(κ), which is the length scale on which the surface wave attenuates parallel to the substrate plane. Given (4.1), the impedance boundary condition (2.12) yields the dispersion relation (κ 2 − Ω 2 ) 1/2 = Ω 2 g/A. With this dispersion relation, the decay constraint and (2.7) imply Ω <Ω.
Our use of the homogenized model is consistent only as long as it predicts surface wavelengths long compared to the spacing. To simplify the discussion below, we assume A 1, whereby the latter condition is κ 1. The dispersion relation then gives g 1/ , and so v , i.e. δ 2 . In that regime, r v and the dispersion relation can be written with the shape factor σ negligible in the expression (2.8) forΩ.
In particular, the homogenized model can be used to study surface waves under the infinitemetasurface critical-coupling condition, since then v 1/2 . For that case, figure 2a depicts the real and imaginary parts of κ as a function of frequency, assuming for later comparison a wave in the x direction:î =ê x and κ = κ x . Consistently with the condition stated above for the homogenization approximation to hold, the dispersion relation 'folds' and κ is small for all frequencies. As κ → 0, one branch joins with the sound line κ = 1/2 Ω, whereas on the other branch Ω →Ω. In the limit (κ, Ω) → (0,Ω), we have b → ∞ and the limiting state propagates to infinity, i.e. it is not a surface wave. In fact, this limiting state is nothing but the solution already found in §3b describing perfect absorption of a normally incident plane wave. Analysis of (4.2) in the critically coupled case v 1/2 reveals two distinguished regimes. The first, defined by the scalings κ − 1/2Ω 5/6 and Ω −Ω 1/3 , corresponds to the dispersion curve separating from the sound line as the resonance frequency is approached. In that regime, p 1/ and b 1/ 2/3 , thus the waves propagate over large distances compared to the freespace wavelength but also decay away from the substrate on a scale larger than the free-space wavelength. In the second regime, defined by the scalings κ 1/2 and Ω −Ω 1/2 , viscosity enters the dominant balance and turns the dispersion curve around. In that regime, p , b 1/ 1/2 , thus the surface waves propagate and decay on length scales comparable to the free-space wavelength. For low-loss metasurfaces such that v is at most of order , i.e. v = O( ), the homogenization model fails as it predicts surface waves with κ 1. Similarly, the homogenization model is not  valid in scattering problems, whether involving a finite or infinite metasurface, where such shortwavelength surface waves are locally excited. Instead, we must return to the multiple-scattering formulation, where the quasi-periodicity condition becomes p n ∝ exp{iκ · r n }. In particular, in the electronic supplementary material (see footnote 1), we show that in the limit → 0, with A 1, the dispersion relation for an inviscid metasurface is well approximated by whereΩ 2 I = K + K 2 σ is the inviscid resonance frequency and S 1 (κ) and S 2 (κ) are the absolutely convergent lattice sums in which Λ = {(n 1 , n 2 ) ∈ Z | A 1/2 (n 1êx + n 2êy )} is the infinite square lattice of neck positions, Λ * = {(n 1 , n 2 ) ∈ Z | (2π/A 1/2 )(n 1êx + n 2êy )} is the corresponding reciprocal lattice and the dash in (4.4a,b) says to omit the zeroth lattice vector. In contrast to the critical-coupling regime, the dispersion relation in the inviscid case is real valued and anisotropic. In figure 2, it is plotted for κ traversing the boundary of the reduced Brillouin zone [29]. We see that the inviscid and criticalcoupling cases, with all other parameters equal, agree for small κ and close to the sound line.
In the inviscid case, however, the dispersion curve does not fold but rather reaches κ 1, where the dispersion curve is extremely flat, indicating small group velocity. In that flat-band regime, the surface wavelength is comparable to the subwavelength periodicity of the metasurface; the decay of the surface wave away from the substrate is also on that length scale, i.e. b 1; and, from (4.3), the dispersion relation can be approximated explicitly (S 1 (κ) + S 2 (κ)). (4.5) In the inviscid case, the propagation length p is infinite. Analysing the propagation length of lowloss metasurfaces, v = O( ), is challenging as it entails deriving a complex dispersion relation starting from the multiple-scattering theory. It is clear from the above, however, that in that case

Large metasurfaces (a) Tuned metasurface, plane-wave forcing
In this section, we present results for finite metasurfaces based on numerical solutions of the multiple-scattering model, and also of the homogenized model where relevant. We first consider the scenario where a plane wave is normally incident on a finite 'tuned' metasurface. By tuned, we mean that the loss parameter δ is tuned to the critical-coupling condition derived in §3b assuming a plane wave normally incident on an infinite metasurface. Figure 3a shows the absorption efficiency D/N as a function of the dimensionless frequency Ω, for the indicated values of N. Values calculated based on the multiple-scattering model are seen to be in good agreement with BEM solutions of the homogenized model. Also shown is the theoretical approximation for an infinite metasurface, calculated using (2.14), (3.2) and (3.3). By design, the efficiency of the infinite metasurface is unity at resonance, corresponding to perfect absorption. The finite metasurfaces are more efficient but, clearly, limited by their size. Qualitatively, their efficiency profiles approach that for an infinite metasurface slowly, nonuniformly (slower below the resonance frequencyΩ) and non-monotonically. Furthermore, in contrast to the infinite-metasurface case, for large metasurfaces the cavity-pressure distribution can be appreciably non-uniform. Figure 3b shows distributions calculated using the multiplescattering model for √ N = 45, at resonance, alongside distributions obtained from BEM solutions of the homogenized model and using (3.4) to define continuous cavity-pressure distributions based on the macroscale pressure field.
Consider now the convergence with increasing metasurface size √ N of the efficiency of finite tuned metasurface to their infinite counterparts. Figure 3f shows the relative difference in efficiency between the finite and infinite cases as a function of √ N, for two values of Ω. The lower frequency is Ω = 0.67. We have seen in §4 that around this frequency, an infinite tuned metasurface supports weakly bound, long-wavelength and attenuating surface waves, with the inverse surface-wavelength 2π/κ, normal decay length b and propagation length p all comparable to the order 1/ √ free-space wavelength. In our infinite-metasurface theory, we have not included these surface-wave solutions as they are not bounded. For a large metasurface, these surface waves are excited, approximately, at the periphery and propagate inwards. This picture is consistent with the observation that the relative difference in efficiency remains order unity up to √ N ≈ 20, roughly p /2 for surface waves at that frequency. Between that size and √ N ≈ 40, roughly p , the relative difference decays by an order of magnitude, presumably owing to the surface waves affecting an increasingly smaller proportion of the metasurface. For larger metasurfaces, we see a much slower convergence, roughly like 1/ √ N. The latter scaling is consistent with the effect of surface-wave excitation being limited to a neighbourhood of the periphery. At the higher frequency, Ω = 0.74, where an infinite tuned metasurface does not support surface waves, we observe a generally smaller relative difference in efficiency.

(b) Detuned (low-loss) metasurface, plane-wave forcing
We next consider the scenario where a plane wave is normally incident on a finite 'detuned' metasurface. By detuned, we mean that the loss parameter δ is not tuned to the critical coupling condition for an infinite system. In particular, we shall specifically consider low-loss detuned metasurfaces for which δ is tuned to the single-resonator critical-coupling condition derived in §3a. According to the infinite-metasurface theory, such low-loss metasurfaces are extremely poor absorbers. In that theory, however, surface waves are excluded. In the detuned case, excitation of surface waves at the periphery can give rise to pronounced finite-size effects even for very large metasurfaces. Based on the inviscid multiple-scattering calculations in §4, for |Ω −Ω| , with Ω <Ω, we expect weakly bound surface waves similar to those in the tuned case. For |Ω −Ω| = O( ), we expect surface waves that are strongly bound, short wavelength and hardly attenuate, with surface wavelength 2π/κ 1, normal decay length b 1 and propagation length p 1/ . Note that the homogenization model is not valid when such short-wavelength surface waves are excited. In this regime, surface waves can traverse even a large metasurface multiple times, reflecting and refracting at the periphery. For certain combinations of frequency and metasurface size, the surface waves may form standing-wave patterns over the square metasurface domain, similar to cavity modes of bulk waves. These constitute leaky states that are damped only by refraction at the periphery and the weak material loss.  Figure 4a shows the averaged dissipation per resonator, or absorption efficiency, plotted as a function of frequency, for tuned and detuned metasurfaces of finite size √ N = 25 and their theoretical infinite counterparts. The results for the finite metasurfaces are obtained by numerically solving the multiple-scattering model, while those for the infinite metasurfaces are calculated using (2.14), (3.2) and (3.3). At most frequencies, the absorption efficiency of the detuned metasurfaces is smaller by orders of magnitude relative to the tuned metasurfaces, as one would expect from the infinite-metasurface theory. There is a narrow frequency interval below the resonance frequencyΩ, however, where the detuned metasurfaces exhibit enhanced absorption, including a sequence of sharp peaks where efficiency is comparable and for some peaks even higher than that of the tuned metasurfaces. We attribute these sharp peaks to collective meta-resonances of the metasurface as a whole, corresponding to excitation of the standing-wave leaky states anticipated above. This picture is qualitatively confirmed by figure 4b, which shows cavity-pressure distributions computed at the first nine spectral peaks. These show standing-wave patterns of shorter and shorter surface wavelength as the cut-off frequency is approached, in accordance with the surface-wave dispersion curve shown in figure 2b for an inviscid metasurface.
For the detuned metasurfaces, figure 4c shows the relative difference in absorption efficiency between the finite and infinite cases as a function of √ N, for three values of Ω. The frequencies Ω = 0.73 and Ω = 0.74 are slightly below the resonance frequencyΩ and the nearby surfacewave cut-off frequency. In contrast to the tuned case, we see no sign of convergence of the efficiency of finite detuned metasurfaces with increasing N, at least up to the values of N we can access numerically. The relative difference remains very large (also because of the very small denominator) and oscillates wildly, more so at Ω = 0.74. Clearly, the apparent lack of convergence is a result of the long propagation length of surface waves at these frequencies, while the wild oscillations are associated with the excitation of meta-resonances, here accessed at fixed frequency by increasing size. The third frequency is Ω = 0.75; it is slightly aboveΩ and the nearby surface-wave cut-off frequency. At that frequency, we observe monotonic convergence, the relative difference in efficiency vanishing like 1/N 17/20 and, after √ N ≈ 20, like 1/N 1/2 . The latter scaling is consistent with local diffraction along the periphery (not triggering surface waves). The origin of the former scaling is not clear.

(c) Tuned and detuned metasurface, monopole-source forcing
Consider now the scenario where instead of a normally incident plane wave there is a localized monopole source positioned at unit height above the centre of the metasurface. Here, surface waves can be excited by the ambient field associated with the source acting on the bulk of the metasurface, as well as by diffraction of that ambient field at the periphery of the metsurface. Thus, in this scenario, surface waves are excited also for an infinite metasurface. We note that, even when short-wavelength surface waves are not excited, the homogenized model fails close to the source owing to its height above the substrate being comparable to the spacing. Figure 5a shows the total dissipation as a function of √ N, for tuned metasurfaces at Ω = 0.67 (below cut-off) and Ω = 0.74 (above cut-off) and for a detuned metasurface at Ω = 0.74 (below cut-off). For large tuned metasurfaces, dissipation at the lower frequency, where weakly bound and attenuating surface-waves are excited, is roughly 2.6 higher than at the higher frequency, where surface waves are not excited. Large detuned metasurfaces, supposed to be poor absorbers, absorb more than large tuned metasurfaces, by an order of magnitude at meta-resonance peaks. Furthermore, the total dissipation oscillates wildly with increasing metasurface size. As in the plane-wave scenario, there is no indication of convergence up to the values of N we compute. Figure 5b contrasts the local dissipation near the source for a tuned metasurface with the distributed dissipation over a detuned metasurface excited at a meta-resonance.

Concluding remarks
We have explored the effects of finite size on the absorption characteristics of large Helmholtztype metasurfaces. The most significant effects are a consequence of the excitation of surface waves, and, for low enough material loss, meta-resonances associated with the surface waves forming standing-wave patterns over the metasurface domain. In particular, finite, even very large, low-loss metasurfaces exhibit absorption characteristics which are nothing like their theoretical infinite counterparts. Thus, while infinite low-loss metasurfaces are poor absorbers, large low-loss metasurfaces are efficient absorbers at near meta-resonance frequencies; this is in comparison with similar 'tuned' metasurfaces satisfying the critical coupling condition for an infinite metasurface. The meta-resonances of low-loss metasurfaces are manifested in a series of sharp spectral peaks accumulating near the cut-off frequency for surface-wave propagation, which is close to the resonance frequency of an isolated resonator. This accumulation is understood to occur owing to the quantization of short-wavelength surface waves, analogously to the accumulation of cavity modes at high frequencies. Thus, low-loss metasurfaces enable precise absorption of wave energy at a series of frequencies. This attribute may be exploited for filtering and sensing applications, assuming that the meta-resonance frequencies could be controlled at will. The present study has focused on modelling and numerical exploration. More sophisticated mathematical analysis is necessary in order to precisely predict the frequencies and quality factors of meta-resonances, as well as address several other open problems highlighted by our examples. These include studying the convergence, with increasing metasurface size, of absorption efficiency and other metasurface properties. As we have seen, the rate of convergence can be extremely non-uniform in frequency, and strikingly different for detuned (low-loss) and tuned (critically coupled) metasurfaces. A first step towards an improved understanding could involve asymptotic analysis of the homogenized problem. One relevant limit process is that of a metasurface comparable in size to the free-space wavelength, considered at a frequency and loss levels such that short-wavelength surface waves are excited. It may be easier to first consider a circular metasurface patch; asymptotic approximations could then be extracted either from an exact Wiener-Hopf solution, or, more insightfully, from a WKB approximation of the surface waves matched with local Wiener-Hopf solutions near the metasurface periphery [30], where the surface waves reflect and refract. Such an approach is only valid, however, in cases where the surface wavelengths are much larger than the periodicity of the metasurface; this must be kept in mind in order to avoid misleading predictions. Accordingly, a second step would be to address scenarios where periodicity-scale surface waves are excited. This could be done by a similar approach to that described above, with the significant technical complication of employing a two-scale WKB approximation [31] to describe the surface waves, matched with local discrete-Wiener-Hopf [32] solutions near the periphery. Of course, many generalizations are possible, including varying the shape of the metasurface domain, metasurfaces formed of non-identical resonators, interaction between multiple metasurface patches etc.
Lastly, since surface waves are intrinsic to metasurfaces formed of subwavelength resonators, we anticipate anomalous finite-size effects similar to what we have found for acoustic Helmholtztype metasurfaces in other acoustic, as well as elastic and photonic, metasurface absorbers. To better understand these effects, our explicit, first-principles and computationally straightforward model could be employed in two ways. First, to design simple experiments to probe these finitesize effects in acoustics. Second, as a toy mathematical model for a generic metasurface absorber, ignoring the extensive geometric and physical details encoded in the lumped parameters.
Data accessibility. The data are provided in electronic supplementary material [33]. Authors' contributions. O.S.: conceptualization, formal analysis, writing-original draft, writing-review and editing; R.B.: conceptualization, formal analysis, writing-review and editing. Both authors gave final approval for publication and agreed to be held accountable for the work performed therein.