# An asymptotic hyperbolic–elliptic model for flexural-seismic metasurfaces

## Abstract

We consider a periodic array of resonators, formed from Euler–Bernoulli beams, attached to the surface of an elastic half-space. Earlier studies of such systems have concentrated on compressional resonators. In this paper, we consider the effect of the flexural motion of the resonators, adapting a recently established asymptotic methodology that leads to an explicit scalar hyperbolic equation governing the propagation of Rayleigh-like waves. Compared with classical approaches, the asymptotic model yields a significantly simpler dispersion relation, with closed-form solutions, shown to be accurate for surface wave-speeds close to that of the Rayleigh wave. Special attention is devoted to the effect of various junction conditions joining the beams to the elastic half-space which arise from considering flexural motion and are not present for the case of purely compressional resonators. Such effects are shown to provide significant and interesting features and, in particular, the choice of junction conditions dramatically changes the distribution and sizes of stop bands. Given that flexural vibrations in thin beams are excited more readily than compressional modes and the ability to model elastic surface waves using the scalar wave equation (i.e. waves on a membrane), the paper provides new pathways towards novel experimental set-ups for elastic metasurfaces.

### 1. Introduction

In recent years, there has been an explosion of interest, both scholarly and popular, in the use of so-called *locally resonant metamaterials* in the control of wave propagation. The concept has its origins in optics [1], but has since found application across a broad range of fields [2], including acoustics [3–6], elasticity [7–9], flexural waves in thin plates [10–12] and seismic waves [13,14]. Despite the term *metamaterial* first appearing at the turn of the millennium [15], the development of this relatively new field continues apace. The essential concept is that, using arrays of local resonators, metamaterials can be designed to have properties that do not occur in materials found in nature; the canonical example is negative refraction [16].

At their inception, metamaterials were designed for the control of waves in bulk media. This idea can be extended to surface waves, creating the so-called *metasurfaces* [17]; these metasurfaces have the substantial advantage over metamaterials in that they are relatively easy to manufacture [18]. The majority of work on metasurfaces has focused on electromagnetism [18] and, for example, spoof-surface-plasmons [19]. At present, however, the use of metasurfaces has begun to gain traction in other fields, such as acoustics [20,21] and, very recently, elastodynamics [22,23] and fluid–solid interactions [24]. In comparison to, say, electromagnetism and acoustics, the treatment of metasurfaces in the framework of elastodynamics is far more challenging and offers a greater array of unique phenomena.

Despite the singular challenges associated with mechanical metamaterials, they have a substantial advantage in that the range of scales over which they can be used is large compared with electromagnetic metamaterials: from the order of millimetres [25] for cloaking of small-scale structures to several metres [13] for seismic metamaterials. The second notable feature is the existence of two bulk waves (shear and pressure), which propagate at two different speeds; this phenomena is unique to elastic waves and has no analogue in electromagnetics or acoustics. The existence of two distinct classes of wave adds both challenges and novel features. Elastic media also support surface waves, again travelling at two different speeds, that are localized to the surface of the medium, but do not decay in the directions parallel to the surface. This is in contrast to electromagnetism where surface plasmons [26] decay in all directions. The closest analogue of ‘true’ surface waves in electromagnetism are the so-called spoof-surface-plasmons [19], which are supported on structured surfaces.

Recently, there have been several analytical [22,23], numerical [14,27] and experimental [28] studies on the use of sub-wavelength resonators for the control of surface waves on elastic substrates. Notable features of such resonant arrays include filtering, focusing, rainbow trapping and mode-conversion of surface waves to bulk waves [29].

Interestingly, a recent numerical and experimental study [30] has shown that the hybridization of the two zeroth-order modes in thin plates can be controlled by using symmetric and asymmetric distributions of resonators. The experimental studies [28,30] are incumbent on prior completion of numerical or analytical studies. The numerical approaches are computationally intensive and require extensive GPU computing facilities [14], while complete analytical studies present their own significant challenges [23].

In the present paper, we develop a novel asymptotic formulation, initially introduced in [31], to analyse the interaction of surface waves with a resonant array resting on an elastic half-space. This asymptotic approach is used to analyse the dynamic effects near the surface of semi-infinite elastic bodies and consists of an explicit hyperbolic equation governing the Rayleigh wave propagation as well as a pseudo-static elliptic equation governing decay into the interior. This approach has proven a powerful tool for a variety of near-surface elastodynamic problems, including three-dimensional moving load problems [32] and the effect of adding pre-stress to the elastic half-space [33]. For more general works, the reader is referred [34–36]. Recently, this hyperbolic–elliptic model was applied to a resonant array of rods attached to a half-space [37]; these rods were assumed to undergo only longitudinal vibrations as was the case in previous analytical [23] and numerical [22] studies. However, recent experimental data obtained for a dense forest of pine trees, acting as subwavelength resonators, suggests that the flexural resonances as well as the junction conditions [38] may play an important role; yet, little attention has thus far been given to the junction condition joining the resonators to the substrate.

Here we extend this approach to include the effects of flexural motion in the thin rods and investigate the effect of these flexural resonances as well as the junction conditions coupling the resonators to the elastic half-space. In order to better see the effect of the flexural motion, and because in the frequency ranges considered flexural waves dominate, the longitudinal motion of the resonators is neglected. These new interactions bring new features, including the ability to control the width of the deep sub-wavelength stop bands. This is an important feature as, usually, sub-wavelength resonant arrays are only capable of producing very narrow band gaps [7] and, as a result, it is challenging to make use of these band gaps in the design of filters and other wave control devices. While a non-dimensionalization of this problem would be useful to consider, this manuscript is following from an existing body of work with strong practical motivations; in the literature of such problems, it is more typical to see the results presented in dimensional form.

As is usual for the treatment of such low amplitude bending waves, the effect of the weight of the beams on vertical force into the half-space and the bending motion of the resonators is considered to be negligible compared to the bending forces.

The present paper is be organized as follows. The problem formulation is given in §2 with the formulation for Rayleigh waves in §2b. The sections which follow then each detail the treatment of a different junction condition between the resonators and the half-space. Section 3 considers the resonators as simply supported, with the resonators able to freely rotate but conserve displacements and horizontal stresses. In §4, the resonators are considered as ideally attached at the base to rails which can move freely but matches the bending moment and the gradient at the surface. For the junction considered in §5, we consider the resonators as ideally attached to the half-space and so all variables between the resonator and half-space are fully matched. In each section, the solution obtained from the asymptotic formulation is compared with an implicit evaluation of the full unimodal solution as used in the previous treatment of similar systems [23].

### 2. Statement of the problem

Before considering the effects of the various interface conditions, we begin with a statement of the problem. We consider the propagation of surface waves on an infinite elastic half-space upon which rests an array of thin Euler–Bernoulli beams, as illustrated in figure 1. These surface waves propagate with constant amplitude along the free surface of the half-space and decay, exponentially, into the bulk. We begin by formulating the three-dimensional problem.

#### (a) An array of beams embedded into an elastic half-space

We consider the linearly elastic motion of a three-dimensional half-space − ∞ < *x*_{1} < ∞, − ∞ < *x*_{2} < ∞, 0 ≤ *x*_{3} < ∞. For the usual displacement vector ** u** = [

*u*

_{1},

*u*

_{2},

*u*

_{3}]

^{T}, we consider only motion in the

*x*

_{1}- and

*x*

_{3}-directions and introduce two displacement potentials

*ϕ*and

*ψ*such that

*ρ*is the volume mass density and where

*λ*and

*μ*are Lamé's first and second parameters, respectively, the equations of motion can be expressed as

*t*is time and ${c}_{1}=\sqrt{(\lambda +2\mu )/\rho}$ and ${c}_{2}=\sqrt{\mu /\rho}$ are the longitudinal and shear wave speed, respectively.

The flexural motion of the thin elastic beams that pattern the surface of the elastic half-space and form the resonant array obey the Euler–Bernoulli beam equation

**(**

*U**t*,

*x*

_{3}) = [

*U*

_{1},

*U*

_{2}]

^{T}is the displacement vector parallel to the surface of the half-space

*x*

_{3}= 0, for simplicity

*E*

_{R}is introduced to represent the resonator stiffness

*E*and

*ν*, respectively, are the Young modulus and Poisson ratio of the resonators,

*M*

_{R}is the mass per unit length of the resonators, and

*I*

_{R}is the area moment of inertia of the beams, which for resonators with cylindrical cross-section is given by

*h*is the diameter of a beam. The ends of the beams at

*x*

_{3}= −

*L*are assumed to be free such that the moments and transverse forces vanish

*x*

_{3}= 0), connecting it to the half-space, undefined. These junction conditions will be defined and studied in subsequent sections. In general, the interface between one-dimensional and three-dimensional bodies results in a singularly perturbed problem which requires careful consideration of the boundary layer (e.g. [39]).

For sufficiently thin resonators, the resistance to bending motion is, typically, much smaller than the resistance to compression and, therefore, in the linear regime, the flexural and compressional deformations are decoupled at leading order (e.g. [40]). Moreover, Colquitt *et al*. [23] suggest that, for the case of resonant arrays on thin plates, the effects of flexural and compressional resonances can be considered independently. Resonant arrays composed of purely compressional resonators atop a half-space have previously been considered in [22,23,37]. However, a detailed analysis of the influence of flexural resonances on the propagation of surface waves is lacking. Therefore, in the present paper, we focus on the flexural motion of the resonators and consider pure bending only such that vertical motion on the surface of the half-space does not couple to compressional deformation of the resonators. It is however emphasized that vertical tractions may arise at the surface of the half-space as a result of rotations in the flexural resonators [40], and these are taken into account.

#### (b) Asymptotic formulation for Rayleigh waves

Considering only waves propagating in the positive *x*_{1}-direction leads to a simplification of the governing equations of elasticity and subsequently a simplification to the asymptotic model, for which the full three-dimensional model is given in [32].

Using the displacement potentials (2.1) and governing equations (2.2) from before, we introduce the usual stress tensor ** σ**, where in the case of a free surface, stress free junction conditions are imposed yielding

*ω*

*A*

_{ϕ},

*A*

_{ψ},

*k*and are constants and

*α*and

*β*are given by

*c*

_{R}, can be obtained. The dispersion equation for such waves is linear,

*ω*=

*c*

_{R}

*k*, and the Rayleigh decay constants are subsequently given by

In the parameter regime of interest, where the characteristic wavelength of the travelling wave is sufficiently large compared with both the cross-section and separation of the resonators, then, as in previous treatments of such problems [37,41], the surface stress associated with the resonant array may be distributed in the form

*H*and

*V*are the amplitudes of the horizontal and vertical forces, respectively, and

*l*is the separation between the resonators.

Each type of junction condition treated in this paper results from the matching of variables between the resonators and the surface of the half-space at *x*_{3} = 0. The same approach will be used in the treatment of the moments arising from the resonators. However, we defer this discussion until the need arises in §4.

We now develop the asymptotic model, initially introduced by Kaplunov *et al.* for a near-Rayleigh surface waves (see [32] and references therein for further details). This asymptotic model is a leading-order perturbation of the equations of motion (2.2) around the Rayleigh wave eigensolution above (2.11). A perturbation in the form of a small surface stress is imposed which causes a small deviation from the Rayleigh speed. This expansion is valid in the regime where

*ω*/

*k*∼

*c*

_{R}. From inspection of the surface conditions, this is equivalent to when both of the surface stresses are small. The equations for the bulk then reduce to

*σ*

_{33}= 0, then it is convenient to work with the shear wave potential and the junction condition along the surface

*x*

_{3}= 0 is

*σ*

_{31}= 0, then the junction condition at

*x*

_{3}= 0 is written in terms of the pressure potential such that

*z*= 0, the displacement potentials are related by

It is straightforward to recover the three-dimensional field from the two-dimensional asymptotic model. In particular, it is necessary to introduce a second shear potential *ψ*_{2}, corresponding to the stress *σ*_{32}, and then simply replace any second-order derivatives with respect to *x*_{1} with the two-dimensional Laplace operator, Δ. For instance, the vertical surface relation (2.19) in the full three-dimensional form becomes

The problem under consideration concerns the effect of a resonant array which, as described earlier, can be treated as a distributed stress on a travelling Rayleigh wave; this makes the asymptotic model introduced here an ideal tool to investigate the problem. For completeness, however, the solution to the full problem (cf. §2a) will also be computed numerically and used as a comparison to the asymptotic solution as verification. We shall now proceed to consider the effect of various interface conditions on the propagation of surface waves over the half-space.

### 3. Simply supported beams

We begin by considering the resonators as simply supported at the surface of the half-space. In this system, there is no need to conserve bending angle and there is no overall bending moment at the coupling. Instead only the horizontal displacement and the transverse force require matching, as illustrated in figure 2. We emphasize that, as discussed in §2a, we consider the influence of the resonators under pure bending only and, since there are no rotations at the junction point, the vertical component of traction does not couple into the beam. Note also that, since the beam rotations are not conserved in the half-plane, any rigid body rotations of the beams are also neglected. The junction conditions at the end of the beam atop the half-space are then

*U*

_{1}is the horizontal displacement of the beam and

*H*(

*x*

_{1},

*y*,

*t*) is the horizontal force at the surface of the half-space caused by the resonators.

By solving the beam equation (2.3) with these junction conditions and the free end conditions (2.6) in the usual way, we obtain the relation for horizontal force

*l*is the distance between two adjacent resonators.

We first consider the full unimodal formulation for the system. From (2.8), the surface stresses can be expressed in terms of displacement potentials such that

*U*

_{1}(0) =

*k*(

*iϕ*+

*βψ*) and, by substituting from (3.3), we obtain the full unimodal dispersion equation

For the asymptotic model, we use equation (2.17) along the surface

*r*=

*ω*

^{2}

*c*

^{−2}

_{2}

*k*

^{−2}. For 0 <

*r*≪1, the leading-order term of (3.12) is

*r*, the full unimodal dispersion relation, (3.8) can be expanded in the form

*β*=

*β*

_{R}.

It is clear that the asymptotic dispersion relation (3.11) is much simpler than the full unimodal equation (3.8) obtained previously. Importantly, the asymptotic dispersion relation can be solved explicitly for *k* as a function of *ω* and, therefore, manipulated and interpreted with greater ease. We compare the asymptotic and full unimodal dispersion relations in figure 3 using the parameter values in table 1.

symbol | definition | value |
---|---|---|

l |
lattice spacing | 2 m |

ρ |
half-space density | 13 000 kg m^{−3} |

μ |
half-space shear modulus | 325 MPa |

λ |
half-space first Lamé parameter | 702 MPa |

L |
resonator length | 14 m |

h |
resonator diameter | 0.3 m |

ρ_{R} |
resonator density | 450 kg m^{−3} |

E_{R} |
resonator stiffness | 1.70 GPa |

The homogenization scheme from (2.13) assumes that the wavelength is on a larger scale than the separation of the resonators, or equivalently *kl*≪2*π*. Subtituting *l* from the system parameters in table 1, this condition is *k*≪*π* *m*^{−1} for each of the illustrative figures, figures 3–6. In addition, while the asymptotic method gives a solution which is valid for all *k*, we note that in the full unimodal system *β* will become purely imaginary when *ω*/*k* > *c*_{2} and so no propagating solution can exist in this region of the plot.

Figure 3 shows good agreement between the dispersion curves for the asymptotic and full unimodal approaches; in particular near the Rayleigh sound-line and for lower frequencies, as expected. In addition, the asymptotes created by the flexural resonances of the beams extends the range of validity of the asymptotic approach beyond the regime of small *k* and small stresses where one would expect the method to remain accurate. The flexural resonances associated with these asymptotes create fine deep sub-wavelength band gaps which have been shown to have a range of physical applications (e.g. [10,22,23,28]).

### 4. Beams on a rail

Suppose now that the resonators are held at *x*_{3} = 0 to a rail parallel to the surface of the half-space, which allows the base of the beam to move freely but changes the bending angle of the beam depending on the gradient of the surface. This allows the beam to impart a bending moment on the half-space but no horizontal forcing. There is therefore no need to conserve displacements and the horizontal stress must vanish. The reader is reminded that, as discussed in §2a, we consider the influence of the resonators under pure bending only. However, in contrast to the case of simple support beams considered in §3, the resonators exert a vertical force on the half-space through rotations in the beams. These rotations can be linked to the gradient of the vertical displacement through the kinematic continuity conditions imposed at the junction points. While this junction condition may initially appear somewhat counterintuitive, it can be understood in terms of the so-called *gyroscopic hinges* [42].

To construct the junction conditions for the base of the beam, we must first determine the bending moment *M* in the elastic half-space caused by the flexural deformation of the beam. For a beam with cross-sectional area *A*, the bending moment is

*V*(

*x*

_{1},

*x*

_{2},

*t*) is the vertical force exerted on the half-space by the beams and

*u*

_{3,1}is the gradient of the vertical displacement along the surface. In particular, solving these in the usual way along with the free-end conditions (2.6) and the beam equation (2.3) gives the vertical force

*H*(

*x*

_{1},

*x*

_{2},

*t*), we will use the force distribution (2.13) and define $\hat{V}$ such that

*u*

_{3,1}in terms of the potentials (2.1) and making use of (2.8) the surface stresses on the half-space can be written in the form

To apply the asymptotic model to this system, where there is only a vertical stress, we will make use of (2.19) along the surface of the half-space, whence we obtain

*ψ*

*R*(

*r*)

*R*(

*r*) for

*r*≪1, that is for wave speeds close to that of a pure Rayleigh wave, to obtain

*α*=

*α*

_{R}and

*β*=

*β*

_{R}, i.e. the decay rates are equal to those of pure Rayleigh waves.

As with the simply supported case, the dispersion relation derived from the asymptotic approach is much simpler than the one produced by the complete analysis. Indeed, the asymptotic dispersion equation is again quadratic and therefore yields an explicit expression for *k* in terms of *ω*. A plot of the solutions of the asymptotic dispersion equation along with those of the full unimodal solutions is shown in figure 4; the parameter values used are again those detailed in table 1. Figure 4 shows that, once again, there is good agreement between the asymptotic model and the full unimodal solution, particularly near the Rayleigh line. For all of the modes shown in figure 4, the asymptotic approach predicts the asymptote, the overall shape of the dispersion relation, and the existence of a large band gap for each mode. Indeed, in addition to being a much more elegant approach, it is clear that the asymptotic model is capable of reproducing the main features of the full unimodal solution. Moreover, compared with figure 3, figure 4 exhibits significantly larger band gaps.

### 5. Fully matched beams

We now consider the case where both displacements and rotations (surface gradients), in addition to, horizontal and vertical tractions are coupled at the junction points. Although, as discussed in §2a, we only consider flexural motion of the resonators and neglect the compressional deformation, the vertical stresses in the half-space couple to the flexural modes in the beams through rotations at the junction point, as in §4; this is fully taken into account here. Predictably, this complex coupling leads to a more involved solution for both the beam motion and the surface wave dispersion equation. For simplicity, we also assume an ideal contact between the beams and half-space, neglecting any local contact strain effects.

In this case, the equations of motion for the beams (2.3) together with the free end conditions (2.6) remain unchanged. The junction conditions at the base of the beam *x*_{3} = 0 are a combination of those previously considered (3.1) and (4.3) with each of the quantities illustrated in figure 2 conserved, giving

In contrast to the simply supported and rail conditions discussed in §§4 and 3, we should not expect the full unimodal solution to intersect with the asymptotic solution at the Rayleigh line. In order for the asymptotic and full unimodal solutions to intersect at the Rayleigh line, the beam must behave as though the coupled end with the half-space is free; this imposes requirements on the displacement and rotation at the base of the beam that are incompatible with the required displacement and gradient for the half-space. We emphasize that this feature does not affect the full unimodal solution but does affect the accuracy of the asymptotic solution since the basic assumption underlying the asymptotic model, that stresses tend to zero at the Rayleigh line, is no longer valid.

To formulate the full unimodal problem, we begin by writing the surface stresses in terms of potentials (2.8), and in terms of the distributed load (2.13), to obtain

The asymptotic model of Kaplunov *et al.* [31] consists of two relations which each employ a single surface stress to produce a wave. However, in the present problem, the system has two stresses; one of which is normal to, and the other of which is tangential to, the surface. Therefore, the approach of Kaplunov *et al.* [31] must be adapted to account for the composition of two stresses into a single wave. In particular, we will assume that there will be two waves, each generated by one of the two stress states. Since, for the asymptotic model both waves have a speed which is a small perturbation from the Rayleigh speed, it can also be assumed that these two waves travel at the same speed. Each of these two waves will have a corresponding pair of displacement potentials, *ψ*_{1} and *ϕ*_{1}, and *ψ*_{2} and *ϕ*_{2} such that by linear superposition the overall displacement potentials *ϕ* and *ψ* are given by a linear combination such that *ϕ* = *ϕ*_{1} + *ϕ*_{2} and *ψ* = *ψ*_{1} + *ψ*_{2}. Each pair of potentials then corresponds to one of the surface relations from (2.17) and (4.11)

As with the previous cases, we note that the asymptotic dispersion equation (5.9) is quadratic in *k* and can therefore be solved explicitly. While more difficult to interpret than the simply supported asymptotic solution, it is still a substantial improvement on the full unimodal dispersion relation. Figure 5 shows the dispersion curves for corresponding to both the asymptotic (5.9) and full unimodal (5.5) solutions. We observe that the asymptotic model predicts the main features of the dispersion curves, including prediction of the existence and general shape of the solution branch for each quasi-periodic mode and also giving a close approximation of the intersections with the Rayleigh line. The key feature of the plot, and one present in both the asymptotic and full unimodal results, is the band gaps. This system is notable for the large range of frequencies which do not have a corresponding real wavenumber solution and will, therefore, not propagate through the system. While it is difficult to locate where the band gaps are from the full unimodal dispersion relation, they can be approximated by the intersections between the asymptotic model and the shear wave line. It is emphasized that, once again, this feature is absent in the previous analyses [22,23,37] as the flexural resonances of the beams were neglected.

This system also highlights some of the limitations of this model. In particular, and unlike for the simple supported and rail junctions conditions that involve a single stress, the asymptotics and the full unimodal solution are not guaranteed to intersect at the Rayleigh line. This is due to the way the stresses are treated in the two methods. In the asymptotic model, there is no interaction between the two stresses and there will be Rayleigh line intersections only at

*α*and

*β*with

*α*

_{R}and

*β*

_{R}and use the surface relation (2.20) to obtain

*α*=

*α*

_{R}and

*β*=

*β*

_{R}and substituting these into the expansion (5.13) leads to the dispersion relation

*k*. However, even for simple stresses, the dispersion equation arising from the aforementioned stress interactions will, in general, be a quartic polynomial.

Figure 6 shows the dispersion curves associated with the full unimodal solution (5.5) and those arising from the asymptotic approach accounting for the stress interaction (5.9); this should be compared with figure 5. Figure 6 illustrates that accounting for the stress interaction forces the model to coincide with the full unimodal dispersion relation at the Rayleigh line. That said, comparing figures 5 and 6 we observe that the dispersion curves arising from the asymptotic approach where the stress interaction is taken into account does not show a significant improvement for the rest of the plot, away from the Rayleigh line, over the standard asymptotic model. This indicates that, for this system, the stress interaction does not have a significant effect on the overall behaviour of the waves away from the Rayleigh line.

Again, comparing figure 3 with figures 5 and 6, it is apparent that the fully matched junction condition allows for the creation of significantly larger deep sub-wavelength band gaps compared with the simply supported junction condition. Therefore, the effect of the flexural resonances on the overall behaviour of the system can be increased by an appropriate choice of junction condition.

A subtlety associated with applying the asymptotic model to multi-stress systems is that the asymptotic approach initial developed in [31] requires that the overall stress is small. It follows that for multi-stress systems, all individual stresses must be small in for the approach to be applicable. If one stress is always large while the other is small then the asymptotic model cannot be valid, even if the result crosses the Rayleigh line.

Finally, away from the Rayleigh line, the dispersion curves approach asymptotes associated with the flexural resonances of the beams; this is a known feature of sub-wavelength resonant arrays (e.g. [11]). For the previous two cases of simply supported and rail junction conditions, the asymptotic model correctly reproduces the asymptotes. This led to the asymptotic model remaining accurate to a good degree even away from the Rayleigh line. By contrast, the fully matched junction (both uni- and multi-stress) does not accurately reproduce the necessary asymptotes. In particular, the asymptotes of the full unimodel solution are solutions of the transcendental equation

*k*becomes large. Therefore, while this is a powerful technique for when both stresses are small, the results cannot be relied upon away from this condition.

Figure 6 also shows the mode shapes associated with each branch of the dispersion curves; the colour scheme shows the magnitude of the total displacement. Animated versions of these figures can be found in the electronic supplementary material. The modes were computed, with the finite element software COMSOL Multiphysics^{®}, using the, using the exact formulation for an infinite periodic array of resonators as described in §2a. It can be observed that each mode can be associated with the fundamental modes of the flexural resonators; this is particularly apparent from the animated mode shapes that can be found in the electronic supplementary material. Indeed, it is clear that the flexural motion of the beams dominates the behaviour of the system in this regime and emphasizes that the flexural resonances can play an important part in such problems where slender beams rest atop an elastic half-space.

Moreover, figure 7 shows the vertical components of the displacement in the beam (*U*_{3}) and half-space (*u*_{3}). It is clear that, in terms of compressional motion, the beams move as a rigid body in the vertical direction and, therefore, the dynamic effects of the compressional motion of the beams do not materially contribute to the overall behaviour of the system in this regime. This observation is particularly evident in the animated versions of figure 7, which are available in the electronic supplementary material.

### 6. Concluding remarks

An array of flexural resonators attached to the surface of an elastic half-space is analysed using an explicit model for the Rayleigh wave. The paper generalizes previous considerations using both full unimodal [23] and asymptotic solutions for an elastic half-space in the case of an array of compressional resonators [37].

Three types of junction conditions are considered, and for each the asymptotic solution is verified using full unimodal solutions for the half-space. The first junction condition to be considered treats the resonators as simply supported at the surface of the half-space. Matching of the beam equation yields a horizontal stress only. The resulting dispersion curves exhibit a series of very fine band gaps typical of deep sub-wavelength arrays; the narrow nature of these band gaps means that making use of these stop bands in practical applications such as filtering and mode-conversion is extremely challenging. For example, deep sub-wavelength arrays cannot be used to design broadband filters.

The second considered junction condition treats the resonators as being supported by a freely moving rail at the surface of the half-space. In contrast to the previous system, this yields a vertical stress only. The effect of this junction condition is to widen the band gaps in the spectrum. This ability to increase the width of the stop bands through a judicious choice of the junction conditions is an important feature that would allow, for example, the implementation of broadband filters in the deep sub-wavelength regime. We also note that this effect will increase the contribution of the dynamic flexural behaviour of the resonators on the overall response of the system. While sub-wavelength resonant arrays have found extensive use in the literature (e.g. [7]), little attention has thus far been given to the junction condition joining the resonators to the substrate. Indeed, one application of the present paper is in the analysis of the interaction of surfaces waves with dense pine forests and recent experimental data suggest that the flexural resonances as well as the junction conditions may be an important factor [38].

For both the simply supported and rail junction conditions, the asymptotic formulation shows good agreement with the full unimodal solution. The asymptotic solution obtained predicts both a coincidence with the full unimodal solution at the Rayleigh line, and also resonances arising from the solution to the beam equation. Moreover, the asymptotic solution has the distinct advantage of being much simpler and easier to manipulate than the full unimodal solution.

The final junction condition considered in this paper fully matches all variables at the junction. Unlike the previous two cases, this results in both a vertical stress and a horizontal stress, requiring an alteration to the existing asymptotic formulation to account for both stresses While this asymptotic solution closely matches with the full unimodal solution, it does not accurately predict the coincidence with the Rayleigh line or the location of the system resonances. A further addition to the asymptotic model, which accounts for the coupling of stresses at the junction point, produces a solution which does accurately predict coincidence with the full unimodal solution with the Rayleigh line, but still does not accurately predict resonances.

All of the systems considered produce band gaps in the region around the beam resonances, with the largest band gaps produced by the beam on a rail and fully matched junction conditions.

This is an important feature as, usually, sub-wavelength resonant arrays are only capable of producing very narrow band gaps [7] and, as a result, it is challenging to make use of these band gaps in the design of filters and other control devices. However, as we have demonstrated, it is possible to control the width of these deep sub-wavelength band gaps by adjusting the junction condition joining the resonators to the half-space. Indeed, the beam on a rail and full matched junction conditions lend themselves to feasible designs of mechanical structure that exhibit large band gaps and are capable of controlling the propagation of surface waves. In particular, the beam on a rail junction condition is associated with gyroscopic resonators, as discussed in [42], and which can be designed to exhibit the desired properties; the fully match condition is, in some sense, the most natural junction condition to impose as it most closely resembles the conditions of ideal contact.

This demonstrates the potential for such systems to be used in controlling the propagation of surface waves.

The explicit asymptotic formulation has the potential to be applied in more sophisticated systems where a full unimodal solution may be difficult or impossible to obtain since it reduces the vector problem in full linear elasticity to a scalar problem along the surface. It has been shown that for the systems considered the asymptotic formulation easily produces simple explicit solutions which match closely with the full unimodal solution.

The formulation discussed involves a surface scalar problem identical to the equation for transverse forces applied to an elastic membrane. Currently, experimental set-ups to model elastic half-spaces can be large and expensive [29] and this formulation gives scope for simpler experiments on membranes which model the behaviour at the surface of an elastic half-space.

This system could also be further developed by considering the near-surface effect of coupled stresses. We also remark that the effect of microstructure, including non-local effects, have previously been considered for near-surface wave motion [43] and, as might be expected, are asymptotically secondary for treating Rayleigh-type waves.

### Data accessibility

This work does not include any experimental data. All computational results presented as examples in the paper are reproducible using the analytical methods detailed in the article.

### Author's contributions

All the authors contributed equally to both, the writing of this paper, and the results presented herein.

### Competing interests

The authors are aware of no competing interests.

### Funding

D.J.C. gratefully acknowledges the financial support of the EPSRC (UK) through programme grant no. EP/L024926/1. P.T.W also acknowledges Keele University for supporting their PhD studies.

### References

- 1.
Shelby RA, Smith DR, Schultz S . 2001 Experimental verification of a negative index of refraction.**Science**, 77–79. (10.1126/science.1058847) Crossref, PubMed, Web of Science, Google Scholar**292** - 2.
Kadic M, Bückmann T, Schittny R, Wegener M . 2013 Metamaterials beyond electromagnetism.**Rep. Prog. Phys.**,**76**126501 . (10.1088/0034-4885/76/12/126501) Crossref, PubMed, Web of Science, Google Scholar - 3.
Liu Z, Zhang X, Mao Y, Zhu Y, Yang Z, Chan CT, Sheng P . 2000 Locally resonant sonic materials.**Science**, 1734–1736. (10.1126/science.289.5485.1734) Crossref, PubMed, Web of Science, Google Scholar**289** - 4.
Li J, Chan CT . 2004 Double-negative acoustic metamaterial.**Phys. Rev. E**,**70**055602 . (10.1103/PhysRevE.70.055602) Crossref, Web of Science, Google Scholar - 5.
Movchan A, Guenneau S . 2004 Split-ring resonators and localized modes.**Phys. Rev. B**,**70**125116 . (10.1103/PhysRevB.70.125116) Crossref, Web of Science, Google Scholar - 6.
Guenneau S, Movchan A, Pétursson G, Ramakrishna SA . 2007 Acoustic metamaterials for sound focusing and confinement.**New J. Phys.**, 399. (10.1088/1367-2630/9/11/399) Crossref, Web of Science, Google Scholar**9** - 7.
Zhou X, Liu X, Hu G . 2012 Elastic metamaterials with local resonances: an overview.**Theor. Appl. Mech. Lett.**,**2**041001 . (10.1063/2.1204101) Crossref, Google Scholar - 8.
Milton GW, Willis JR . 2007 On modifications of Newton's second law and linear continuum elastodynamics.**Proc. R. Soc. A**, 855–880. (10.1098/rspa.2006.1795) Link, Google Scholar**463** - 9.
Achaoui Y, Khelif A, Benchabane S, Robert L, Laude V . 2011 Experimental observation of locally-resonant and Bragg band gaps for surface guided waves in a phononic crystal of pillars.**Phys. Rev. B**,**83**104201 . (10.1103/PhysRevB.83.104201) Crossref, Web of Science, Google Scholar - 10.
Colombi A, Roux P, Rupin M . 2014 Sub-wavelength energy trapping of elastic waves in a metamaterial.**J. Acoust. Soc. Am.**, EL192–EL198. (10.1121/1.4890942) Crossref, PubMed, Web of Science, Google Scholar**136** - 11.
Williams EG, Roux P, Rupin M, Kuperman W . 2015 Theory of multiresonant metamaterials for*A*_{0}lamb waves.**Phys. Rev. B**,**91**104307 . (10.1103/PhysRevB.91.104307) Crossref, Web of Science, Google Scholar - 12.
Xiao Y, Wen J, Wen X . 2012 Flexural wave band gaps in locally resonant thin plates with periodically attached spring–mass resonators.**J. Phys. D: Appl. Phys.**,**45**195401 . (10.1088/0022-3727/45/19/195401) Crossref, Web of Science, Google Scholar - 13.
Brûlé S, Javelaud E, Enoch S, Guenneau S . 2014 Experiments on seismic metamaterials: molding surface waves.**Phys. Rev. Lett.**,**112**133901 . (10.1103/PhysRevLett.112.133901) Crossref, PubMed, Web of Science, Google Scholar - 14.
Colombi A, Roux P, Guenneau S, Gueguen P, Craster RV . 2016 Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances.**Sci. Rep.**, 19238. (10.1038/srep19238) Crossref, PubMed, Web of Science, Google Scholar**6** - 15.
Walser RM . 2001 Electromagnetic metamaterials. In*Proc. SPIE 4467, Complex Mediums II: Beyond Linear Isotropic Dielectrics, San Diego, CA*, 9 July 2001. Bellingham, WA: SPIE. Google Scholar - 16.
Veselago VG . 1968 The electrodynamics of substances with simultaneously negative values of*ε*and*μ*.**Sov. Phys. Usp.**, 509–514. (10.1070/PU1968v010n04ABEH003699) Crossref, Web of Science, Google Scholar**10** - 17.
Maradudin AA . 2011**Structured surfaces as optical metamaterials**. Cambridge, UK: Cambridge University Press. Crossref, Google Scholar - 18.
Chen HT, Taylor AJ, Yu N . 2016 A review of metasurfaces: physics and applications.**Rep. Prog. Phys.**,**79**076401 . (10.1088/0034-4885/79/7/076401) Crossref, PubMed, Web of Science, Google Scholar - 19.
Pendry J, Martin-Moreno L, Garcia-Vidal F . 2004 Mimicking surface plasmons with structured surfaces.**Science**, 847–848. (10.1126/science.1098999) Crossref, PubMed, Web of Science, Google Scholar**305** - 20.
Ni X, Wu Y, Chen ZG, Zheng LY, Xu YL, Nayar P, Liu XP, Lu MH, Chen YF . 2014 Acoustic rainbow trapping by coiling up space.**Sci. Rep.**, 7038. (10.1038/srep07038) Crossref, PubMed, Web of Science, Google Scholar**4** - 21.
Zhou C, Yuan B, Cheng Y, Liu X . 2016 Precise rainbow trapping for low-frequency acoustic waves with micro Mie resonance-based structures.**Appl. Phys. Lett.**,**108**063501 . (10.1063/1.4941664) Crossref, Web of Science, Google Scholar - 22.
Colombi A, Colquitt D, Roux P, Guenneau S, Craster RV . 2016 A seismic metamaterial: the resonant metawedge.**Sci. Rep.**, 27717. (10.1038/srep27717) Crossref, PubMed, Web of Science, Google Scholar**6** - 23.
Colquitt D, Colombi A, Craster R, Roux P, Guenneau S . 2017 Seismic metasurfaces: sub-wavelength resonators and Rayleigh wave interaction.**J. Mech. Phys. Solids**, 379–393. (10.1016/j.jmps.2016.12.004) Crossref, Web of Science, Google Scholar**99** - 24.
Skelton E, Craster R, Colombi A, Colquitt D . 2018 The multi-physics metawedge: graded arrays on fluid-loaded elastic plates and the mechanical analogues of rainbow trapping and mode conversion.**New J. Phys.**,**20**053017 . (10.1088/1367-2630/aabecf) Crossref, Web of Science, Google Scholar - 25.
Bückmann T, Thiel M, Kadic M, Schittny R, Wegener M . 2014 An elasto-mechanical unfeelability cloak made of pentamode metamaterials.**Nat. Commun.**, 4130. (10.1038/ncomms5130) Crossref, PubMed, Web of Science, Google Scholar**5** - 26.
Ritchie R . 1957 Plasma losses by fast electrons in thin films.**Phys. Rev.**, 874–881. (10.1103/PhysRev.106.874) Crossref, Web of Science, Google Scholar**106** - 27.
Palermo A, Krödel S, Marzani A, Daraio C . 2016 Engineered metabarrier as shield from seismic surface waves.**Sci. Rep.**, 39356. (10.1038/srep39356) Crossref, PubMed, Web of Science, Google Scholar**6** - 28.
Colombi A *et al.*2017 Enhanced sensing and conversion of ultrasonic Rayleigh waves by elastic metasurfaces.**Sci. Rep.**, 6750. (10.1038/s41598-017-07151-6) Crossref, PubMed, Web of Science, Google Scholar**7** - 29.
Colombi A, Craster RV, Colquitt D, Achaoui Y, Guenneau S, Roux P, Rupin M . 2017b Elastic wave control beyond band-gaps: shaping the flow of waves in plates and half-spaces with subwavelength resonant rods.**Front. Mech. Eng.**, 10. (10.3389/fmech.2017.00010) Crossref, Google Scholar**3** - 30.
Rupin M, Roux P, Lerosey G, Lemoult F . 2015 Symmetry issues in the hybridization of multi-mode waves with resonators: an example with Lamb waves metamaterial.**Sci. Rep.**, 13714. (10.1038/srep13714) Crossref, PubMed, Web of Science, Google Scholar**5** - 31.
Kaplunov YD, Kossovich LY . 2004 Asymptotic model of Rayleigh waves in the far-field zone in an elastic half-plane.**Dokl. Phys.**, 234–236. (10.1134/1.1753618) Crossref, Web of Science, Google Scholar**49** - 32.
Kaplunov J, Prikazchikov DA . 2017 Asymptotic theory for Rayleigh and Rayleigh-type waves.**Adv. Appl. Mech.**, 1–106. (10.1016/bs.aams.2017.01.001) Crossref, Web of Science, Google Scholar**50** - 33.
Khajiyeva L, Prikazchikov D, Prikazchikova L . 2018 Hyperbolic-elliptic model for surface wave in a pre-stressed incompressible elastic half-space.**Mech. Res. Commun.**, 49–53. (10.1016/j.mechrescom.2018.07.006) Crossref, Web of Science, Google Scholar**92** - 34.
Kaplunov J, Zakharov A, Prikazchikov D . 2006 Explicit models for elastic and piezoelastic surface waves.**IMA J. Appl. Math.**, 768–782. (10.1093/imamat/hxl012) Crossref, Web of Science, Google Scholar**71** - 35.
Ege N, Erbas B, Chorozoglou A, Kaplunov J, Prikazchikov DA . 2017 On surface wave fields arising in soil-structure interaction problems.**Procedia Eng.**, 2366–2371. (10.1016/j.proeng.2017.09.253) Crossref, Google Scholar**199** - 36.
Nobili A, Prikazchikov DA . 2018 Explicit formulation for the Rayleigh wave field induced by surface stresses in an orthorhombic half-plane.**Eur. J. Mech. A/Solids**, 86–94. (10.1016/j.euromechsol.2018.01.012) Crossref, Web of Science, Google Scholar**70** - 37.
Ege N, Erbaş B, Kaplunov J, Wootton P . 2018 Approximate analysis of surface wave-structure interaction.**J. Mech. Mater. Struct.**, 297–309. (10.2140/jomms) Crossref, Web of Science, Google Scholar**13** - 38.
Roux P *et al.*2018 Toward seismic metamaterials: the METAFORET project.**Seismol. Res. Lett.**, 582–593. (10.1785/0220170196) Crossref, Web of Science, Google Scholar**89** - 39.
Kozlov V, Maz'ya VG, Movchan AB . 1999**Asymptotic analysis of fields in multi-structures**. Oxford, UK: Oxford University Press on Demand. Google Scholar - 40.
Kaplunov JD, Kossovitch LY, Nolde E . 1998**Dynamics of thin walled elastic bodies**. New York, NY: Academic Press. Google Scholar - 41.
Slepyan L . 1967 The strain wave in a bar with vibration-isolated masses.**Mech. Solids**, 57–64. Google Scholar**2** - 42.
Nieves M, Carta G, Jones I, Movchan A, Movchan N . 2018 Vibrations and elastic waves in chiral multi-structures.**J. Mech. Phys. Solids**, 387–408. (10.1016/j.jmps.2018.07.020) Crossref, Web of Science, Google Scholar**121** - 43.
Chebakov R, Kaplunov J, Rogerson G . 2016 Refined boundary conditions on the free surface of an elastic half-space taking into account non-local effects.**Proc. R. Soc. A**,**472**20150800 . (10.1098/rspa.2015.0800) Link, Google Scholar