To facilitate the interpretation of results presented in this article, we first recall the particularity of Rayleigh waves in terms of ellipticity and we present the specific analysis of the ground motion polarization, thanks to data recorded at 10 m from the impact in Fig. 4. Finally, we propose an approximate plate model that qualitatively reproduces the observed effects.

### Ellipticity of Rayleigh waves

In the left part of Fig. 4, we present the velocity data (x, y, z) recorded at 10 m from the impact. In the right part, we show the related plots of particle motion in horizontal (x–y) plane and vertical (x–z) plane. In this case, the signal is strongly polarized in the horizontal plane and more exactly in the x-direction i.e. perpendicularly to the side of the mesh. This information is essential because it corroborates that most of the energy of the source is converted into energetic surface waves. This ellipticity may complicate the interpretation of results and the effectiveness of the grid of holes. As for velocity dispersion, the Rayleigh wave motion is frequency dependent and such an ellipticity function can be represented by the Fourier transform ratio of the horizontal and the vertical components of motion versus frequency. The shape of this function could be “flat”, if the velocity of compressional and shear waves for first superficial strata do not vary brutally versus depth at each interface. The preliminary seismic tests carried out on the virgin ground, confirm the flatness of this function in the range of 2 to 11 Hz, in accordance with the geological data.

### Numerical results

Let us now describe some numerical simulations, which provide some phenomenological understanding of wave phenomena at stake in the experiments led by the Ménard company. Much has been said about control of light, sound, water, or shear (SH) waves31,32,33,34,35,36 using the rich behavior encapsulated by the dispersion curves of bi-periodic structures, modeled by a Helmholtz equation

$Δu+k2u=0$

, with

$Δ=∂2∂x2+∂2∂y2+∂2∂z2$

the Laplacian, up to minor changes in the normalization of material parameters, and choice of boundary conditions (e.g. Dirichlet or Neumann for clamped or freely vibrating inclusions in the context of SH waves). However, when one considers elastic waves, governed by the vector Navier equations, which for an isotropic homogeneous medium take the following form in the frequency domain (1):

$(λ+2µ)∇(∇·U)−µ∇×(∇×U)=−ρω2U$

(1)

with ω the angular wave frequency (in unit of rad.s−1), λ and µ the Lamé parameters (Pa), and ρ the density of the elastic medium (kg.m−3), things become more complicated. Indeed, the unknown U(x,y,z) is a 3-component displacement field, and

$∇=(∂∂x,∂∂y,∂∂z)$

denotes the gradient with partial derivatives in all three coordinates, whereas . and × denote the usual scalar and cross products in Cartesian coordinates. When the elastic medium is homogeneous, one can still decouple the Navier equation(1) in Helmholtz-type equations, but if one considers stress-free inclusions in the elastic medium, it is no longer possible to reduce the analysis to scalar partial differential equations (PDEs), as shear and pressure waves do couple at inclusions’ boundaries. This makes difficult correspondences between electromagnetic and elastodynamic metamaterials.

There is nevertheless, the simplified framework of the Kirchhoff-Love plate theory37 that allows for bending moments and transverse shear forces to be taken into account via a fourth-order PDE for the out-of-plane plate displacement field. This plate theory is a natural extension of the Helmholtz equation to a generic model for flexural wave propagation through any spatially varying elastic medium thin compared with the typical wave wavelength. This PDE is deduced from an asymptotic analysis in (1), that amounts to taking a limit of small plate thickness in out-of-plane coordinate z, which leads to37,38

where D = Eh3/(12(1 − ν2)) is the so-called plate flexural rigidity, h its thickness, E the Young’s modulus and ν the Poisson ratio, which are such that E = λ(1 + ν)(1 − 2ν)/ν. Moreover, u is the third component of U.

Kirchhoff-Love equation(2) offers a very convenient mathematical model to draw analogies between the physics of electromagnetic and mechanical metamaterials. One can then transfer earlier knowledge gained in photonic or platonic crystals to large scale seismic metamaterials. However, while the Helmholtz equation can, with appropriate notational and linguistic changes, hold for acoustic, electromagnetic, water or out-of-plane elastic waves and so encompasses many possible applications, the Kirchhoff-Love plate theory is dedicated to the analysis of flexural waves (it does not model propagation of in-plane elastic waves in platonic crystals). Nonetheless, (2) can be used as a simplified, but meaningful, model for Rayleigh waves in structured soils10. Correspondences between models for scalar waves in photonic crystals and surface elastic waves in structured soils are unveiled when one recasts Eq. (2) as per:

$(∂2∂x2+∂2∂y2+Ω)(∂2∂x2+∂2∂y2−Ω)u=0$

(3)

where

$Ω2=12(1−ν2)ρω2Eh2$

Upon inspection, (3) looks like a factor of two Helmholtz equations with spectral parameters k2 = Ω and (ik)2 = Ω with i2 = −1. Physically this means that one has coexistence of propagating and evanescent flexural waves in the plate.

To model the soil structured with boreholes, we consider a plate with same geometric and soil parameters as in the experiments and with thickness 5 m, which contains an array of 23 stress-free inclusions. This is an approximate model where the moment and the transverse force are both equal to zero. This means we set the following two boundary conditions on each air hole (written in polar coordinates for simplicity as we consider circular holes):

$D∂2u+∂r2+Dv(1∂ur∂r+∂2ur2∂θ2)=0$

(4)

$D∂∂r(∂2u∂r2+1∂ur∂r+∂2ur2∂θ2)+1rD(1−v)∂∂θ(1r∂2u∂r∂θ−∂2ur2∂θ)=0$

(5)

where (4) enforces vanishing moment and (5) enforces vanishing transverse force.

It is possible to show using high-frequency homogenization21 that the solution u satisfies the following effective dispersion relation

$Ω∼Ω0−Tij/(2Ω0)κiκj$

where Tij is a rank-2 tensor encompassing the (frequency dependent) effective elastic parameters, κi is a component of the Bloch wavenumber and Ω0 is the standing wave frequency at the Brillouin zone edge 0, −π/2, π/2 depending on the location in the Brillouin zone about which the asymptotic expansion originates. For symmetry reasons, only diagonal components Tii
are non- zero. Besides from that, when T11T22 < 0 the effective dispersion relation describes some hyperbolic type wave propagation i.e. waves propagate like in the geometrical ray optics theory19, and this happens notably at a frequency around 8 Hz, see Fig. 5a and b, where a point source gives rise to an image inside and a faded one outside a flat slab with effective parameters T11 = 31 and T22 = −14 whereas T11 = T22 = 1 outside the slab, in accordance with same simulation for the array of air holes, see Fig. 5c and d for plots of Re(u) and |u|, respectively. Note that elastic parameters of soil are embedded in Ω for which a Young modulus E = 0.153 GPa, a Poisson ratio ν = 0.3 and a density ρ = 1800 kg m−3 are considered, as well as a plate thickness h = 5 m. This theoretical prediction for a hyperbolic-type metamaterial is in accordance with Fig. 3d, where both the red and green curves display a local maximum around 8 Hz, which coincides with a minimum for blue curve. On the other hand, when |T11| << |T22|, one observes some cross like wave propagation, which happens at frequency around 3 Hz, and is akin to the endoscope effect discussed in19. We have checked that these theoretical predictions still hold for the full elasticity case by solving the full Navier equations (1) for a plate with thickness 15 m, with stress free boundary conditions on cylindrical inclusions of diameter 2 m and height 15 m (considering cylindrical inclusions of height 5 m does not alter the plots at 8 Hz), and all other boundaries of the domain. Unwanted reflections on all lateral sides of the computational domain are avoided using 3D perfectly matched boundary layers implemented as described in25. One can clearly see the lensing effect on the 3D plot of the modulus of the vector displacement field |U| (Fig. 5e), as well as on 2D plots of a slice of the real part Re(U3) (Fig. 5f) and modulus |U3| (Fig. 5g) of its vertical component taken close to the upper boundary of the plate. Note that we consider in 3D some plate of thickness 15 m to make sure we generate ‘Rayleigh-like’ waves and not Lamb waves anymore like in the 2D flexural wave case, see25 for similar computations for the case of stiff inclusions in soil clamped to a bedrock.

As one can see, this dynamic effective anisotropy is essential in seismic metamaterials as it enables to mould the flow of surface elastic waves almost at will, as was theoretically shown in the case of thin plates19, with an experimental proof of concept for an acoustic source generating hyperbolic-type Lamb wave patterns in a periodically pinned plate at kHz frequencies39. Importantly, homogenized models of thin and thick plates have much in common, and the former provide useful guidance for the latter26. For instance, analogies between models of thin and thick plates makes possible the design of seismic cloaks via artificial anisotropy25,40.