Seismic metasurfaces: Sub-wavelength resonators and Rayleigh wave interaction

We consider the canonical problem of an array of rods, which act as resonators, placed on an elastic substrate; the substrate being either a thin elastic plate or an elastic half-space. In both cases the flexural plate, or Rayleigh surface, waves in the substrate interact with the resonators to create interesting effects such as effective band-gaps for surface waves or filters that transform surface waves into bulk waves; these effects have parallels in the field of optics where such sub-wavelength resonators create metamaterials, and metasurfaces, in the bulk and at the surface respectively. Here we carefully analyse this canonical problem by extracting the dispersion relations analytically thereby examining the influence of both the flexural and compressional resonances on the propagating wave. For an array of resonators atop an elastic half-space we augment the analysis with numerical simulations. Amongst other effects, we demonstrate the striking effect of a dispersion curve that transitions from Rayleigh wave-like to shear wave-like behaviour and the resultant change in displacement from surface to bulk waves.


Introduction
Metamaterials, as synthetic composite materials with a structure such that they exhibit properties not usually found in natural materials, now form a major emerging research area that barely existed before 2000; in fact, the term "Metamaterial" itself was first used in 1999. Since then, the area has grown extensively and shows little sign of slowing down. The key point is that materials can be designed to have, say, a negative refractive index as predicted by Veselago (1968) and later byPendry (2000); and subsequently fabricated by Smith et al. (2000), which is impossible in naturally occurring materials. The first metamaterials were developed in optics and electromagnetism and relied upon having simultaneously negative permittivity and permeability, this was made physically possible using a microstructured periodic medium consisting of sub-wavelength resonators such as split-ring resonators (Lagarkov et al., 1997;Pendry et al., 1999). It has subsequently been realised that these ideas can also be profitably utilised to create acoustic or elastic (negative density and negative shear or bulk modulus) metamaterials (Craster and Guenneau, 2012;Deymier, 2013;Liu et al., 2000;Fang et al., 2006;Kadic et al., 2013). Similarly, although metamaterials were initially developed for bulk media, one can also create microstructured surfaces that act as metasurfaces (Maradudin, 2011) with, again, most of the activity centered around electromagnetic waves and surface plasmons. Given that the array of resonators we introduce here modifies the surface wave properties we choose to call this a metasurface rather than a metamaterial although there is no strict definition, as yet.
Much more recently seismic metamaterials have begun to be considered, although here the challenges are substantial: not only are the Rayleigh waves, which are of primary interest, surface waves, but the underlying system of equations are the full vector equations of elasticity. Nonetheless there have been attempts to modify the local properties of the ground through the addition of inclusions, or resonators, of a different material at a sub-wavelength scale. The different types of inclusions, resonant or non-resonant, determine the properties and the performance of the structured medium or metamaterial/ metasurface. Turning to the non-resonant case first, Brûlé et al. (2014) show with both large-scale experiments and theory that the soil properties can be critically affected by periodic arrangements of boreholes with a spacing of about a metre; in the Bragg scattering regime these induce bandgaps that can be used for seismic protection. More recently, resonant sub-wavelength scatterers placed on top of an elastic substrate have been investigated and we are motivated by recent geophysical experiments (Colombi et al., 2016b) that have demonstrated that natural forest trees can act as a metasurface for frequencies between 30 and 100 Hz. Quite remarkably these experiments show a distinct reduction in the transmission of waves over a broad range of frequencies; these frequencies appear to lie in bandgaps created by local resonances between trees and elastic waves in the substrate. These large-scale experiments can be interpreted through laboratory-scale experiments and theory (Colombi et al., 2014;Williams et al., 2015;Yoritomo et al., 2016) that have, for simplicity, utilised elastic plates, not the full elastic system; notably, these papers neglect the flexural deformations of the resonators. In contrast, in the present paper, we account for both the flexural and compressional deformations of the resonators and demonstrate that the flexural resonances significantly impact the spectrum of the metamaterial in certain regimes. Although we have focused upon the geophysical setting, elastic waves are also very important in ultrasonics and surface acoustic wave devices, and the Bragg scattering aspects of surface arrays as phononic crystals have been explored, i.e., in, say, Achaoui et al. (2011). Locally resonant structures have also been shown to give rise to so-called super-wide pseudo-directional band gaps in platonic crystals (Xiao et al., 2012); such band gaps are formed by the coalescence of Bragg and resonance band gaps.
Given the interest in this emerging area, and the developing applications such as surface to bulk elastic wave filters (Colombi et al., 2016a), there is a need for the fundamental solutions to canonical problems involving periodic arrays of sub-wavelength resonators atop elastic substrates. Here we provide the required analytical and theoretical background for these seismic metasurfaces by considering arrays of these resonators attached to elastic substrates, thereby considering the interaction with elastic Rayleigh waves and for comparison and completeness we also consider the elastic plate. We do so in a two-dimensional setting as the essential concepts are uncluttered by excessive algebra, and the solutions we obtain allow for interpretation and are, on occasion, completely explicit.
The mathematical tools for treating periodic arrays are extensive: linear partial differential operators, quasi-periodic Green's functions, and integral transforms based upon Fourier series and we use a variety of these approaches freely to obtain the cleanest solutions; these theoretical approaches are complimented by illustrative numerical calculations based on the Spectral Element Method. We begin in section 2 by considering the array of resonators attached to an elastic plate, this is complementary to Williams et al. (2015), and we obtain explicit dispersion rela- tions and explore the relative importance of flexural and compressional resonances in the array and how they interact with the plate waves. This naturally leads in the full elastic situation of an array atop an elastic half-space, section 3 investigates this and again explicit results emerge that are used to interpret and predict metamaterial/metasurface phenomena. Finally, we draw together some concluding remarks in section 4.

Thin elastic plates
We begin by considering an infinite periodic array of beams (which support both flexural and compressional waves, whereas rods only support compressional waves) on a thin elastic plate, as shown in figure 1; the boundary of the elementary cell is indicated by the dashed lines. We use Cartesian coordinates x 1 , x 2 orientated along the thin elastic plate, and perpendicular to it with u 1 , u 2 being the associated components of the displacement vector. It is convenient to introduce two sets corresponding to the two different components within each cell. The first defines the flexural foundation G n = {x : |x 1 | < (n + 1)/2, x 2 = 0}, whilst the second corresponds to the resonator R n = {x : x 1 = n, 0 < x 2 < } for n ∈ Z. We note that G n ∩ R m = ∅ and G n ∩R n = {x : x = (n, 0)}.
The equations of motion for time-harmonic waves of radian frequency ω are where β 4 = ρh/(EI), α 2 = ρ/E, E is the Young's modulus of the foundation, I is the area moment of inertia of the foundation, and h is the thickness of the foundation; the subscript R denotes the corresponding properties for the resonators. The amplitudes of the horizontal and vertical forces and moments at the base of the resonators are denoted by F n , V n , and M n respectively; the Dirac delta function, δ, and its derivative, δ , are understood in the distributional sense.
The tops of the resonators are free leading to boundary conditions of the form Imposing continuity of forces and moments at the base of the resonators yields for x ∈ R n ∪ {x : x 2 = 0} where S is the cross-sectional area of the resonators. These boundary conditions are then supplemented with kinematic equations corresponding to continuity of displacements and rotations at the base of the resonator u 1 (n, 0 − ) = −u 1 (n, 0 + ), u 2 (n, 0 − ) = u 2 (n, 0 + ), ∂ ∂x 1 u 2 (n, 0 − ) = − ∂ ∂x 1 u 1 (n, 0 + ), (2.4) where x 2 = 0 − belongs to the closure of G n and x 2 = 0 + belongs to the closure of R n . The negative signs in the first and final continuity conditions are a result of the choice of coordinate system and forcing orientation. Introducing the Fourier transform and its inverse as respectively, together with the Bloch-Floquet quasi-periodicity conditions on the forces and moments (·) n = (·)e iξn where ξ ∈ (−π, π), the fields are expressed in the form where we have restricted ourselves to x 2 = 0 and the second argument is therefore omitted. Making use of the Poisson summation formula (Lighthill, 1958), the fields can be further reduced to infinite sums. Setting x 1 = 0 leads to a 3 × 3 linear system for the displacements and rotations, the solvability condition for which yields the dispersion equation det σ − I = 0, (2.8) where I is the identity matrix and ω 2 being the dispersion equations for a uniform infinite membrane and plate respectively, (2.10e) It now remains to evaluate the infinite sums appearing in the matrix problem.

Evaluation of the infinite sums
We start by considering the simplest summation which, using Poisson summation, can be expressed as (2.11) The integrand has four simple poles located at γ = ±β √ ω and γ = ±iβ √ ω. In order to evaluate the integral we take the usual semi-circular contour closing it in the upper half-plane for n ≥ 0 and the lower half-plane otherwise. For n ≥ 0 the contour is indented such that the pole at γ = −β √ ω lies outside the contour and γ = β √ ω is enclosed by the contour. The converse is true for n < 0. Computing the residues we find that For the remaining summations we find that (2.14) and finally . (2.15)

The dispersion equation
Now that the 3 × 3 matrix σ is expressed in a finite number of terms, it is straightforward to compute the determinant (2.8) and hence obtain the dispersion equation The coefficients A n (ω) are given in A and contain the information about the material and geometrical properties of the resonators. The first parenthesised term in equation (  positive real roots for ω and therefore can be ignored; the second term is the dispersion equation for flexural waves in a uniform thin plate. The remaining parenthesised term in equation (2.16) is a cubic polynomial in cos ξ and therefore has exact closed-form solutions. An example of a typical dispersion diagram is shown in figure 2, the corresponding material and geometrical parameters are detailed in table 1. The dispersion curves for compressional and flexural waves in the foundation, without resonators, are also shown in figure 2. The grey dashed horizontal lines are associated with the flexural resonances of the beams; the corresponding boundary value problem is that of an Euler-Bernoulli beam with one end clamped at x 2 = 0 and the remaining end free at x 2 = . In this case, the resonances satisfy the transcendental equation (Graff, 1975) We note that the density of the resonances reduces with increasing frequency, tending toward the expected constant value (Graff, 1975). The compressional resonances, indicated by the thick solid grey lines in figure 2, are associated with the spectrum of longitudinal waves in a clamped-free thin elastic rod; the natural frequencies are thus ω (r) n = (2n − 1)π/(2 α R ). Comparing figures 2a and 2b we observe that, away from the flexural resonances (indicated by the dashed grey lines), the spectrum appears virtually unchanged by the flexural deformations of the resonators. This effect can be attributed to the contrast in rigidities between the plate and resonator: the plate is roughly 4.5 times as rigid as the resonators, and so the flexural deformations of the resonators couple weakly to the plate. Figure 3 shows the dispersion curves for a system where the base plate is thinner with h = 3 × 10 −3 m, but otherwise identical to that considered earlier; all the other parameters are as detailed in table 1. In this configuration, the rigidity of the plate is approximately half that of the resonators and results in the flexural deformations of the resonators coupling strongly with the flexural waves in the base plate. Indeed, comparing figures 2 and 3, we note that the flexural resonances (indicated by the dashed grey lines) begin to play a more important role. For the thicker more rigid plate, the compressional resonance of the resonators give rise to band gaps. However, for the thinner more flexible plate, the confinement of the dispersion curves between adjacent flexural resonances gives rise to additional band gaps as well as regions of negative group velocity (c.f. the region near ω = 0.225 and ξ = π in figure 3); this negative group velocity can be associated with so-called Double-negative acoustic metamaterials, as discussed by Li and Chan (2004). In previous works, these flexural resonances have been neglected but, as we see here they can play an important role in certain regimes.
This effect can be understood in terms of the forces and moments exerted at the junctions between the plate and resonators. In particular, equations (2.10) give the force (resp. moment) per unit displacement (resp. rotation) exerted by the resonators on the plate. The compressional deformations of the resonators couple to the plate by means of the force V (2.10a), whilst M θ (2.10b), M u (2.10c), F θ (2.10d), and F u (2.10e) couple the flexural deformations to the plate. Examining equations (2.10), we observe that decreasing the thickness of the plate, and hence the area moment of inertia, increases the magnitude of the moments M θ and M u , which couple the flexural deformations of the resonators to the plate.
It has been observed experimentally by Roux et al. (In press) that, for sufficiently thin plates, sharp transmission bands appear in the band gaps created by the longitudinal resonances of the resonators. Figure 4 shows a magnified region inside the first band gap of the spectra illustrated in figures 2 and 3; Panel (a) corresponds to the thicker more rigid plate, whilst panel (b) is associated with the thinner more flexible plate. As can be observed from figure 4, these narrow transmission bands (corresponding to the solid black curves) are associated with the flexural resonances of the resonators. Specifically, the confinement of the dispersion curves between adjacent flexural resonances results in regions of the spectrum with very small, but non-zero, group velocities (the gradients of the dispersion curves); this in turn results in very narrow transmission bands inside the band gaps created by the longitudinal resonances. As we transition from the more rigid plate, c.f. figure 4a, to the more flexible plate, c.f. figure 4b, the group velocity and the band width increase and therefore one would expect the effect to be more noticeable in an experimental setting. However, we emphasise that this effect persists in both the rigid and flexible case and serves to further confirm both the analysis presented here and the earlier experimental results; additionally, this effect highlights that the effects of the flexural resonances can be important, particularly for flexible plates.

Rod-like resonators
From the previous discussion it would appear that, for sufficiently rigid plates, the flexural component of the resonators interacts weakly with the plate and its primary influence is seen as a discrete set of resonances on the dispersion diagram. With this in mind, we now examine the case where the flexural interaction is considered negligible. Under this assumption we may set

An elastic half-plane
We now move to the more geophysically relevant, but involved, problem of a linear array of resonators on an elastic half-plane, as shown in figure 5. In particular, we are interested in the propagation and control of Rayleigh waves, that is, waves that propagate along the surface of the half-space and decay exponentially into the bulk. We begin by formulating the problem as a doubly periodic two-dimensional square array of resonators on an elastic half-space and introduce x 1 x 2 x 3 L Figure  the sets R = {x : x 1 = 0, 0 < x 2 < , x 3 = 0} and H = {x : |x 1 | < L/2, x 2 < 0, |x 3 | < L/2}, which correspond to the resonator and half-space in the elementary cell. It was shown, in §2 as well as in previous experimental (Colombi et al., 2014) and numerical (Williams et al., 2015) investigations, that the contribution from the bending deformations of the resonators does not significantly affect the overall behaviour of the system away from flexural resonances. Therefore, and in order to ease exposition, only the longitudinal motion of the resonators will be considered; the approach used here is equally applicable to the case when the flexural deformation of the resonators is included. With this being the case, the boundary value problem for the displacement amplitude field u = u(x) is expressed as where C is the elastic stiffness tensor of the half-space and n is the outward unit normal. The first equation (3.19a) is the equation of motion for bulk elastic waves in the half-space, whilst (3.19b) is the equation of motion of longitudinal waves in the resonators. Finally, equation (3.19c) represents the balance of forces at the junction between the resonators and half-plane; continuity of u 2 across x 2 = 0 at x 1 = x 3 = 0 is also imposed. For the isotropic half-space considered here, it is convenient to employ the Helmholtz decomposition and express the displacement field in terms of the usual scalar and solenoidal vector potentials (Achenbach, 1984) u(x) = ∇ϕ(x) + ∇ × ψ(x). (3.20) The compressional and shear potentials both satisfy Helmholtz equations with Λ 2 c = ρω 2 /(λ + 2µ), Λ 2 s = ρω 2 /µ, and λ and µ being the Lamé parameters of the elastic half-space. The necessary components of the tractions, corresponding to the boundary conditions (3.19c), are then σ 12 = µ(2φ ,12 − 2ψ ,11 − Λ 2 s ψ), (3.22a) Having formulated the problem in three spatial dimensions, which was necessary in order to ensure correct dimensionality, we now focus on the problem analogous to that considered in §2: a linear array of resonators on a half space. To this end, and motivated by the results of §2, we impose Bloch-Floquet quasi-periodicity and search for solutions in the form of Fourier series with α 2 n = (k − 2πn/L) 2 − Λ 2 c and β 2 n = (k − 2πn/L) 2 − Λ 2 s . The branch cuts are chosen such that (α n ) > 0 and (β n ) > 0 in order to ensure decay as x 2 → −∞. Physically, the above Ansätze correspond to surface waves propagating parallel to the x 1 -axis over the surface of an elastic half-space with a periodic array of slender rods.
Combining (3.19c) and (3.23), we find Multiplying both sides of (3.24) by e −i(k−2πm/L)x1 (m ∈ Z) and integrating over (x 1 , x 3 ) ∈ (−L/2, L/2) 2 yields i.e. the analysis of sub-wavelength control of surface waves, it is sufficient to consider a single mode expansion such that φ m = φ 0 δ m0 and ψ m = ψ 0 δ m0 . It is clear that solutions of this form satisfy (3.25) and, hence, (3.26). In this case, the dispersion equation reduces to where we recognise the left-hand side as the usual Rayleigh dispersion equation (Achenbach, 1984), and where V (ω) = Sω √ EP tan( ω P/E) is the vertical force exerted on the half-plane by the resonators. We have also introduced the normalised variables ξ = k/Λ s and r 2 = Λ 2 c /Λ 2 s = 1/(2 + λ/µ), and 0 < r 2 < 3/4. Although the reduced dispersion equation is transcendental in both ξ and ω and does not, in general, permit closed form solutions, it can be expressed as a sixth order polynomial in ξ 2 allowing the roots to be determined efficiently using the various fast algorithms available for finding polynomial roots.  Figure 6 shows the dispersion curves for a typical configuration; the material parameters are detailed in table 2. The dispersion curves, which are solutions of (3.27), are shown as solid black lines. These curves correspond to combinations of frequency, ω, and wavenumber, ξ, for which a surface waves exists; that is, waves that propagate over the surface of the half-space and decay exponentially into the bulk. These waves arise as a result of the periodicity of the array present on the boundary of the half-space and are distinct from the usual Rayleigh waves that exist on  The geometrical and numerical parameters used to produce the dispersion curves for the half-plane system shown in figure 6. free surfaces; it is natural, therefore, to refer to such waves as Rayleigh-Bloch waves (Porter and Evans, 1999;Colquitt et al., 2015). In addition to the dispersion curves, we also indicate the compressional resonances of the resonators, shown as the horizontal dashed lines in figure 6. The remaining dashed lines of increasing slope are the curves corresponding to Rayleigh waves, bulk shear waves, and bulk compressional waves in the homogeneous elastic half-space respectively; we call these dispersion lines the Rayleigh/shear/compressional wave sound-line by analogy with the terminology used in electromagnetism (see, for example, (Joannopoulos et al., 2011)).

Dispersive properties and critical points
Previous investigations of similar systems were based on purely numerical simulations (Colombi et al., 2016a,b); with the dispersion equation (3.27) in hand, we can now examine several interesting features of the dispersion curves offering additional physical insight. As noted in the earlier papers (Colombi et al., 2014;Williams et al., 2015) on elastic plates and from the numerical simulations of Colombi et al. (2016b) for a half-space, the onset of band gaps coincide with the longitudinal resonances of the resonators. Here we observe that the upper boundary of the band gaps coincide with the intersection of the dispersion curves with the shear sound-line. These points are indicated by the red circles in figure 6 and correspond to the case when ξ = 1 (ω/k = v s ), in which case, equation (3.27) reduces to L 2 √ µρ + √ 1 − r 2 S √ EP tan( ωα R ) = 0 which has the solutions (3.28) The longitudinal resonances of the resonators lie at ω (r) n = (2n − 1)π/(2 α R ). Thus, in this regime where the pass bands are bounded from below by the resonances of the resonators and from above by the intersections of the dispersion curves with the shear sound-line, the band gaps have a constant width defined by (3.29) Since 0 < r 2 < 3/4, the argument of arctan in equation (3.29) is real and positive; hence n > 0 and therefore band gaps always exist for finite parameters. Equation (3.29) provides a convenient formula for tuning the width of the band gap: In order to maximise the band gap, the factor α R and the argument of arctan in (3.29) should be made as small as possible simultaneously. The converse will minimise the width of the stop band. In rough terms, short resonators with high compressional wave-speeds will maximise the stop band width. However, as the material parameters of the resonators also appear in the argument of arctan, some care is required when tuning the system.

Application to filters, metawedges, and effective band-gap widths
Thus far we have studied infinite periodic structures. However, metamaterials, metasurfaces, and phononic structures often find application as filters, lenses, and polarisers. As has been shown in the electromagnetic (Maradudin, 2011), elastodynamic (Colquitt et al., 2011), and thin plate (Dubois et al., 2013) literature, the dispersive properties of an infinite system can be used to design finite structures that posses interesting properties such as negative refraction, flat lenses, filtering, and cloaking. To illustrate the potential of Rayleigh waves in this regard, we consider their interaction with a finite array of resonators. To this end, and given that we are concerned with the control of surface waves on elastic bodies, the intersection of the dispersion curves with the Rayleigh wave sound-line for an homogeneous elastic half-space without resonators (indicated by the blue dot on figure 6) is of interest. For combinations of wavenumber and frequency corresponding to this point of intersection, a Rayleigh wave incident on the array of resonators from the elastic half-space can generate a Rayleigh-Bloch wave, that is a surface wave whose behaviour is affected by the periodicity, in the array. Conversely, at this point, a Rayleigh-Bloch wave leaving the array can couple into a pure Rayleigh wave in the ambient half-space. With the exception of embedded surface waves, which we do not consider here, the homogeneous half-space does not support localised Rayleigh waves corresponding to points above the Rayleigh wave sound-line; in practical terms, this means that the region between the lower red dot and the blue dot on figure 6 is an "effective stop band" for Rayleigh waves incident on the array from a homogeneous half-space.
The point of intersection between the dispersion curves and the Rayleigh sound-line can be obtained from (3.27) by setting ξ = v r /v s , where v r is the Rayleigh wave speed (see (Graff,

1975), among others). The frequency of intersection is then
where γ = v s /v r . The "effective width" of the first band gap is then (3.31) We now illustrate the dispersive properties of the resonant array atop an elastic half-space using the results from time domain spectral element (SEM) simulations. We consider a linearly elastic, isotropic, homogeneous medium with a linear array of resonators, as depicted in figure 7; the material and goemetrical properties of the system are specified in table 2. The simulations are performed with SPECFEM2D 1 (Komatitsch and Vilotte, 1998), a well known parallel code widely used by the seismological community. The mesh is constructed using quadrilateral elements and the commercial software CUBIT 2 . Perfectly matched layers (Komatitsch and Martin, 2007) (PML) are applied on the bottom and vertical boundaries of the computational domain containing the half-space; the remaining boundaries (top surface and resonators) are traction free. The accuracy of this method has been thoroughly tested in a previous study (Colombi et al., 2016b) against experimental measurements.
Thus far we have considered an infinite array of resonators, however for the numerical simulations (Figs. 8 and 9), a large, finite, number is more convenient. We take an array of 30 (3.32) The distance between the source and the resonant array is such that only Rayleigh waves scatter on the resonators (e.g. Figs. 8 and 9); the remaining classes of waves are absorbed by the PML. Once the ramp function has reached its final value, the simulations continue in a stationary state where only a monochromatic signal with constant amplitude propagates. The snapshots in Figs. 8 and 9 have been taken when the simulation is in this regime. Figures 8 and 9 demonstrate the metamaterial properties of the resonant array in various regimes. Here we have chosen not to present the deformation of the resonators in order to highlight the wavefield in the half-space; the full videos associated with these figures are provided in the supplementary materials 3 . Figure 8c shows the transmission of a surface wave through the array at a frequency corresponding to the lower branch of the dispersion curves, figure 8a. We also note the decrease in wavelength inside the array corresponding to the slower effective wavespeed inside the metasurface as compared with the homogeneous half-space. The shorter wavelength is accompanied by an increase in amplitude of the waves inside the array as well as an increased rate of decay into the bulk. In contrast, figure 8b shows the case when the frequency of the incident wave almost coincides with a resonance of the resonators leading to virtually no transmission.
Moving on to figure 9, the two panels (b) and (c) correspond to the second branch of the dispersion curves, as indicated. It is important to note that this upper branch transitions from the Rayleigh sound-line to the shear sound-line, that is, the Rayleigh-Bloch wave in the resonant array preserves the Rayleigh-like surface wave properties at higher frequencies but will evolve into a shear-like wave as the frequency decreases. The hybrid nature of this branch, sitting between the shear and Rayleigh sound-lines, is responsible for the mode conversion of Rayleigh waves into shear waves and used in the design of so-called resonant metawedge systems (Colombi et al., 2016a).
The point of intersection (marked in red) between the hybrid branch and the shear soundline is of particular interest. At this frequency, waves can propagate past the metasurface but they must only do so as shear waves (note the polarization) as is illustrated in figure 9(c) for a frequency very close to this point. In this regime the transition between Rayleigh waves propagating in the free half-space (to the right hand side of the array) and the hybrid mode in the array happens abruptly resulting in an "apparent bandgap" for surface waves. Therefore the wavefield beneath the resonators, represented in figure 9(c), is characterised by shear waves that propagate away from the surface and into the bulk leaving the left hand side of the array almost untouched. As the frequency increases (figure 9(b)) the wavefield reveals more clearly its hybrid nature. The waves propagating under the resonators are grazing the surface but they are not yet exponentially decaying as Rayleigh waves. In this regime waves can propagate in the resonator but the transmission is still weak. The hybrid mode progressively turns into a Rayleigh wave as the frequency increases toward the crossing point (blue dot in figure 6).

Concluding remarks
Metamaterials barely existed fifteen years ago, but now form a major research area: although initially developed for optics, the field of metamaterial research has expanded rapidly and now includes the development of novel materials for applications in acoustics and elasticity. More recently, the applicability of metamaterials to seismology has sparked the interest of geophysicists in the development of novel methods to control surface waves. Given the interest in this emerging area, there is a need to study the properties of the solutions to fundamental canonical problems.
In this paper we have provided the theoretical and analytical framework necessary for the rigorous study of these seismic metasurfaces that adjust the surface wave behaviour. The framework is introduced through the study of two canonical problems associated with the control of mechanical surface waves. Initially, we study the propagation of surface Bloch waves through an array of resonators on a thin elastic plate. In this case, the resonators are thin elastic rods supporting both compressional and flexural waves. Explicit exact solutions are developed and used to examine the behaviour of the system. In particular, we study the coupling between the symmetric and anti-symmetric modes in the substrate and the resonances of the rods. Interestingly it is found that, in the frequency range of interest and for sufficiently rigid plates, the flexural resonances of the resonators couple very weakly into the substrate with the dominant effect coming from the compressional modes of resonators, which open up band gaps associated with the compressional resonances. For more flexible plates, the flexural deformations of the resonators become important.
The plate system is used as the motivation for the far more challenging problem considered in section 3 where we examine the propagation of surface Bloch waves through an array of resonators resting on a fully-elastic half-space. Here we develop closed form expressions for the dispersion equation and wave-fields in the deep sub-wavelength regime of interest. As for the plate problem, the compressional resonances of the resonators create band gaps in the dispersion curves and analytical expressions for their position and width are provided. In particular, it is shown that band gaps will always exist regardless of the material or geometrical parameters of the half-space and resonators. The notion of "effective band gaps" for surface waves is also introduced and discussed.
In section 3, we also examine the scattering problem associated with a half-line of resonators and use the dispersive properties of the array to examine, the filtering effects of the resonant array.
With the formal framework of seismic metasurfaces now established, the time is ripe to exploit these results and explore the possibilities of extending the existing results for electromagnetism and acoustics to seismology. Indeed, the theoretical framework developed in the present paper is already being used in the design and development of the so-called meta-wedge that is capable of mode-converting destructive seismic surface waves into mainly harmless bulk shear waves (Colombi et al., 2016a).