Control of Rayleigh-like waves in thick plate Willis metamaterials

We explore interactions of elastic waves propagating in plates (with soil parameters) structured with concrete pillars buried in the soil. Pillars are 2 m in diameter, 30 m in depth and the plate is 50 m in thickness. We study the frequency range 5 to 10 Hz, for which Rayleigh wave wavelengths are smaller than the plate thickness. This frequency range is compatible with frequency ranges of particular interest in earthquake engineering. It is demonstrated in this paper that two seismic cloaks' configurations allow for an unprecedented flow of elastodynamic energy associated with Rayleigh surface waves. The first cloak design is inspired by some approximation of ideal cloaks' parameters within the framework of thin plate theory. The second, more accomplished but more involved, cloak design is deduced from a geometric transform in the full Navier equations that preserves the symmetry of the elasticity tensor but leads to Willis' equations, well approximated by a homogenization procedure, as corroborated by numerical simulations. The two cloaks's designs are strikingly different, and the superior efficiency of the second type of cloak emphasizes the necessity for rigor in transposition of existing cloaks's designs in thin plates to the geophysics setting. Importantly, we focus our attention on geometric transforms applied to thick plates, which is an intermediate case between thin plates and semi-infinite media, not studied previously. Cloaking efficiency (reduction of the disturbance of the wave wavefront and its amplitude behind an obstacle) and protection (reduction of the wave amplitude within the center of the cloak) are studied for ideal and approximated cloaks' parameters. These results represent a preliminary step towards designs of seismic cloaks for surface Rayleigh waves propagating in sedimentary soils structured with concrete pillars.


I. INTRODUCTION
In a recent theoretical proposal, 1 control of body waves was numerically demonstrated with a spherical shell consisting of an anisotropic heterogeneous elasticity tensor without the minor symmetries (what falls within the framework of so-called Cosserat media), and an heterogeneous (but isotropic) density. In Ref. 1, it was suggested that a simplified version of this elastodynamic cloak could be implemented to protect objects, nuclear waste or even large scale infrastructures buried deep down in the soil. However, achieving asymmetric elasticity tensors for instance via effective medium approaches remains a challenge. Indeed, the Cosserat theory 2 of elasticity (or micropolar elasticity) incorporates a local rotation of points as well as the translation assumed in classical elasticity, and a a Email: andre.diatta@fresnel.fr; b sebastien.guenneau@fresnel.fr couple stress (a torque per unit area) as well as the force stress (force per unit area). In the isotropic Cosserat solid or micropolar continuum, there are six elastic constants, in contrast to the classical elastic solid which is described by two constants. This makes Cosserat cloaks fairly challenging to engineer. In the present paper, we follow a different route towards seismic wave cloaking for surface Rayleigh waves, without resorting to Cosserat media. Before we move to the core material of this paper, we would like to first recall why soils structured at a meter scale might counteract deleterious action of certain types of seismic waves, what might seem at first glance fairly counter-intuitive. More than a million earthquakes are recorded every year, by a worldwide system of earthquake detection stations, some of which are particularly devastating and cause human casualties, such as the earthquake of magnitude-6.2 that struck Italy on the 22nd of August 2016 at 01:36 GMT, 100 km north-east of Rome, not far from L'Aquila, where a similar earthquake struck on the 9th of April 2009. The propagation velocity of the seismic waves depends on density and elasticity of the earth materials (clearly, it is much different in sedimentary soils and rocks). At the scale of an alluvial basin, such as in L'Aquila, seismic effects involve various phenomena, such as wave trapping, resonance of whole basin, propagation in heterogeneous media, and the generation of surface waves at the basin edge. 3 As noted in Ref. 4, wherein meter-scale inertial resonators were introduced to reflect surface and body seismic waves, due to the surface wave velocity in superficial and under-consolidated recent material (less than 100 to 300 m.s −1 ), wavelengths of surface waves induced by natural seismic sources or construction work activities are shorter than those of earthquake generated direct P (primary, i.e. longitudinal compressional) and S (secondary, i.e. transverse shear) waves (considering the 0.1-50 Hz frequency range), from a few meters to a few hundreds of meters. These are of similar length to that of buildings, therefore leading to potential building resonance phenomena in the case of earthquakes such as in L'Aquila. Other sources of soil vibrations one may wish to attenuate include traffic and construction works, 5 and screening methods have been proposed to achieve that. 6,7 Interestingly, back in 1999 Rayleigh wave attenuation was theoretically and experimentally achieved in marble quarry with air holes displaying kHz stop bands 8 and similar filtering effects in a microstructured piezoelectric for MHz till 1 GHz surface waves. [9][10][11][12][13] In our previous work on control of surface Rayleigh waves in soils structured with borehole inclusions, we used so-called elastic stop bands which are frequency intervals for which waves are disallowed to propagate in an infinite periodic array through Bragg interferences (when the wave wavelength is on the order of the array pitch). With boreholes 0.3m in diameter, with a center to center spacing of 1.7m and the soil parameters on this field experiment, 14 this meant we had a stop partial band around 50 Hz (for only one crystallographic direction). Three rows of boreholes were enough to reduce the energy of the signal by about 30 per cent behind the seismic metamaterial. However, if one were to design a shield for the frequency interval 5 to 10 Hz, it would be necessary to scale up the array of boreholes by a factor of five to ten, which requires a lot of free space around the area, one wishes to protect. To overcome this obstacle, we now propose to use columns of concrete instead of boreholes, so as to achieve subwavelength control of Rayleigh surface waves. We show in figure 8 an example of such a structured soil, which was designed by civil engineers of the Menard company.
Importantly, the contrast between the soil and concrete column's parameters makes it possible to achieve uncommon effective material parameters, such as with strong artificial anisotropy and dispersion, and even artificial inertia and viscosity that can be interpreted as rank-3 tensors in an effective Willis equation. This leads to a markedly enhanced control of surface Rayleigh wave trajectories compared to our earlier work. 14 Indeed, comparisons between numerical results demonstrate that one can achieve some cloaking of surface (Rayleigh-like) wave trajectories, where the wave field is almost unperturbed outside the seismic cloak, whereas it simultaneously nearly vanishes in its center. Bearing in mind that we use soil and concrete elastic parameters for the thick plate and pillars, and that we consider wave frequencies lower than 10 hertz, we claim that our results represent a preliminary step towards implementation of seismic cloaks in civil engineering.
In this paper, we would like to present a novel approach for the design of seismic cloaks, which is based on the Willis-type equations, as proposed in the seminal paper by Milton, Briane and Willis, 15 which investigated form invariant governing equations in physics, notably coordinates transforms in Navier equations ensuring a transformed symmetric elasticity tensor. For the sake of simplicity in the numerical computations and physical discussions, we shall restrict ourselves to elastic waves in thick plates, rather than in semi-infinite media. Such a choice is in compliance with the fact that Rayleigh wave amplitudes decay exponentially (it vanishes within two wavelengths) with depth, becoming approximatively zero by the time a depth of two wavelengths has been reached. However, we consider wavelengths small compared to the plate's thickness, hence surface waves are akin to Rayleigh, rather than Lamb, waves. Indeed, as noted in Ref. 14, the elastic plate model already gives some interesting insight in the physics of seismic waves. The frequency of interest in earthquake engineering is governed by the frequencies of structures for the fundamental mode and the first harmonics. Typically, the fundamental period expressed in seconds, for N-story structure could be approximated by N/10. For our numerical implementations, we have decided to focus on the range of 5 to 10 Hz corresponding to low-rise common building with only few storeys.
The plan of the paper is as follows: We first present some derivation of transformed Navier equations with a symmetric elasticity tensor, for thick plates. The obtained Willis equations are then simplified (removal of rank-3 tensors) and the illustrative case of a cloak with heterogeneous anisotropic Willis equations is numerically solved with finite elements implemented in the COMSOL commercial package. We then explain how one can approximate these anisotropic elastic parameters with an effective medium approach. The touchstone of our homogenization approach is that all tensors are symmetric. We finally propose two designs of soils structured with concrete pillars which are numerically validated with finite elements.

II. ELASTODYNAMIC WILLIS MATERIAL, EQUATIONS OF MOTION
The propagation of elastic waves is governed by the Navier equations. Assuming time harmonic exp(−iωt) dependence, with ω as the angular wave frequency and t the time variable, allows us to work directly in the spectral domain. Such dependence is assumed henceforth and suppressed, leading to where ρ is the density of the (possibly heterogeneous isotropic) elastic medium and u = (u r , u θ , u z ) is the three-component displacement field in a cylindrical coordinate basis x = (r, θ, z). Also, C is the rank-four (symmetric) elasticity tensor with components C ijkl (i, j, k, l = r, θ, z) and: stands for the double contraction between tensors. For example, C : ∇u is the 2-tensor with components (C : ∇u) ij = C ijkl ∇u kl . Let us consider a coordinate change x → x , where x = (r , θ , z ) are stretched spherical coordinates. In general, this leads to a transformed equation 15,16 or, in a more compact way where S, D are 3-tensor fields, ∇ is the gradient in transformed coordinates x and u (x ) = u (r , θ , z ) is a transformed displacement in stretched cylindrical coordinates. Note that the transformed stress is generally not symmetric. Note also that in general A is a matrix field. In order to preserve the symmetry of the stress tensor, one assumes that A is a multiple A = ξ∂x /∂x, of the Jacobian matrix ∂x /∂x of the transformation, where ξ is a non-zero scalar, in which case (2), (3) are said to be Willis-type equations. 15,17

III. APPROXIMATED CLOAKING WITH WILLIS MATERIAL WITHOUT 3-TENSORS
Following the proposal of Pendry et al. of an invisibility cloak via transformation optics, 18 we consider a specific geometric transformation, in cylindrical coordinates with θ = θ, z = z. Such a transformation blows the axis {(r, θ, z), r = 0} of a vertical solid cylinder of radius r 2 , into the cylinder {(r, θ, z), r = r 1 }, while confining the whole solid cylinder {(r, θ, z), r ≤ r 2 } into the following hollow cylinder {(r, θ, z), r 1 ≤ r ≤ r 2 } of inner and outer radii r 1 and r 2 respectively. In the sequel, we plug the transformation (4) into the Navier equations to derive the corresponding transformed equations (2), (3). From here on, we drop (neglect) the 3-order tensors S and D (i.e. we neglect inertia and viscosity features in the modified Willis equations) and hence use an approximate version of the Willis-like Equations. Note that in this case, C obviously still has its minor and major symmetries, but is anisotropic and heterogeneous, with coefficients in cylindrical coordinates The stretched density ρ is now a tensor field of order 2, in particular, it is anisotropic and inhomogeneous. Namely, it has non-constant different eigenvalues and is diagonal in cylindrical coordinates with components as follows Although all the (minor and major) symmetries of C are preserved, the anisotropy in the azimutal direction is infinite at the inner boundary of the cloak, as C θθθθ /C r r r r = r r −r 1 4 recedes to infinity as quickly as 1 r −r 1 4 when r approaches r 1 . However C r r r r and C zzzz are of the same order and both tend to zero at the same rate at the inner boundary r =r 1 . Meanwhile, at the same boundary, the off-diagonal components C r r θθ , C θθr r , C rθrθ , C r θθr , C θr r θ , C θr θ , C θθθθ , C θθzz , C zzθθ , C θzθz , C θzzθ , C zθθz , C zθzθ , all become infinite at the rate 1 r −r 1 , whereas C r r zz , C zzr r , C r zr z , C r zzr , C zr r z , C zr zr , C θzθz , C θzzθ , C zθθz , C zθzθ , tend to zero. One remarks in the meantime that, the eigenvalues ρ r r and ρ z z of the density in the r and z directions both vanish, whereas in the azimuthal direction ρ θ θ becomes infinite at the inner boundary of the cloak. The material requirement for vanishing density along the vertical direction ρ z z can be relaxed for Rayleigh waves, since their amplitude goes to zero within two wavelengths (this is fortunate as it means we won't have to structure the soil in the vertical direction when we approximate the ideal cloak parameters through homogenization approach). However, variation of transverse anisotropic density (in the horizontal xy-plane) according to (7) is a strong requirement for a fully operational cloak.
Let us point out at this stage that other approaches to elastodynamic wave cloaking exist, notably for out-of-plane shear 37 and in-plane coupled shear and pressure 38 waves in pre-stressed media, for pressure waves in pentamode solids that behave in a way similar to anisotropic fluids 39 and flexural waves in thin plates. 40 We shall come back to the latter in the sequel as we are particularly interested in plate waves in this paper.

A. Rayleigh-like waves in thick plates
In this paper, we apply the approximate Willis cloaking discussed above, to Rayleigh waves in particular in the framework of seismic metamaterials. Rayleigh waves, unlike Love waves, are surface waves that do not need a layer to guide them. The former surface waves exist in semi-infinite homogeneous isotropic media with a stress -free boundary. In our present work, we consider plates whose thickness is larger than the wavelength of the surface wave, which are akin to Rayleigh waves. When we make the geometric transform in the Navier equations, the stress-free boundary conditions also undergo the mapping. Hence one has to apply the transformation to the Navier Equations coupled with the stress-free boundary condition with n the (outward) normal at the points on the surface z = z 0 . In our case, the horizontal plane z = z 0 stands for the soil-air interface of the thick plate. In the resulting Willis-type equations discussed in the previous sections, the approximated technique is again applied coupled with n · (C : ∇ u) = 0 on the transformed boundary z = z 0 , where n is the normal at the transformed boundary.

B. Numerical implementation of Rayleigh-like waves in thick plates with soil parameters using finite element method
Aiming at applications particularly in civil engineering such as seismic protection, we base the numerical model corresponding to the cloak described above on soil parameters. Namely, we use the Lamé coefficients λ = 87.7 * 10 6 Pa, µ = 58.5 * 10 6 Pa and density ρ = 1800 kg/m 3 . Although the approach developed here is general, we chose to concentrate our numerical tests on the range of 5 to 10 Hz corresponding to low-rise common building with only few storeys.
The cloak consists of a vertical cylindrical hollow region of height 50 m, with an annulus-like basis of inner and outer radii r 1 = 10 m and r 2 = 21.6 m respectively and the z-axis as its main axis. It surrounds a homogeneous isotropic cylindrical domain of radius r 1 (the cloaked region) made of soil. Likewise, the cloak itself is located inside a cylindrical homogeneous isotropic medium (the ambient space) also made of soil and surrounded by a cylindrical shell of inner and outer radii r 3 = 50 m and r 4 = 66.6 m respectively, which is filled with an anisotropic heterogeneous absorptive medium acting as a (reflectionless) perfectly matched layer (PML). The top plane is a stress-free surface, which is the interface between the soil and the air.
The implementation of this model in the finite element package COMSOL MULTIPHYSICS requires the use of Cartesian coordinates. Fortunately, 40 of the 3 4 spatially varying entries of the transformed elasticity tensor vanish identically. We mesh the computational domain using 1,462,046 tetrahedral elements, 54,726 triangular elements, 1,660 edge elements and 41 vertex elements, as shown in Fig. 1(b).
A point force located at (8.3m,36.7m,16.7m), -that is, it lies on the soil-air free surface at a distance of 37.6 meters from the axis of the cloaked region, -oriented along the direction (0,0,1) generates a Rayleigh-like wave at frequency 5.5 Hz Fig(2 (A))-which corresponds to a wave wavelength ∼ 33 m that is greater than half of the basin depth (25 m), so this surface wave still behaves like a plate waveand at frequency 7.9 Hz Fig(2 (B))-which corresponds to a wave wavelength ∼ 23 m strictly lower than half of the basin depth-so it does not feel the bedrock. Following predictions in Ref. 26, one expects that some Rayleigh-like waves generated by specific frequencies of the source will create resonances (trapped modes) within the cloaked region as in Figure 2 (a), wherein a dipole mode can be observed (the inner cylinder wobbles): This occurs at a countable set of eigenfrequencies corresponding to the problem of the Neumann (stress-free) cylindrical cavity in the center of the cloak, as we will investigate in section VI. At a source frequency away from the cavity trapped modes, one can clearly see that we achieve a good seismic protection (Fig. 2 (b)), see also Fig. 7 where a quantitative study is carried out in the frequency range 3 Hz to 10 Hz. Elastic parameters for a cloak with a Willis-like heterogeneous anisotropic transformed medium seem unachievable in practice with current civil engineering techniques. We therefore present in the sequel a homogenization route towards achievable Willis-like seismic cloaks.
At this stage, let us point out that an alternative route to control of Rayleigh waves consists in applying tools of conformal optics as proposed by Ulf Leonhardt, 20 that in the context of transformational elastodynamics lead to spatially varying isotropic elastic parameters. 21 Although this might seem an easier way towards cloaking, this requires structuring soil with pillars softer than soil, which is actually the opposite case to what we do in the sequel.

IV. HOMOGENIZATION APPROACH FOR APPROXIMATE WILLIS MEDIUM WITHOUT RANK-3 AND RANK-2 TENSORS
Our homogenization approach to achieve required rank-4 and rank-2 tensors proceeds in two steps. Firstly, we follow the proposal of Greenleaf and coworkers 22 to apply the two-scale convergence method of G. Allaire and G. Nguetseng 23,24 in the design of approximate multi-layered cloaks, we derive a set of effective parameters for a multi-layered elastic cloak, which is akin to homogenized parameters obtained by the G-convergence approach in Ref. 25, except that these parameters remain dependent upon the macroscopic variable. We consider an alternation of concentric layers of respective elasticity tensors C 1 and C 2 and of respective densities ρ 1 and ρ 2 so that the overall elasticity tensor and density can be written where x is the position vector, and r = x = x 2 1 + x 2 2 is its Euclidean norm and 1 I is the indicator function of the interval I (I is typically the unit cell size along the direction of periodicity). The spatially varying tensor C η is periodic on I = [0,1] and η defines the periodicity (the thinner the layers, the smaller η, the larger the number of layers). The governing equation reads: When η tends to zero, it is shown in Ref. 25 that the solution u η to the above Navier equation two-scale converges towards u 0 solution to the homogenized Navier equation: with an effective (symmetric) elasticity tensor and < f > (x) =< f (x, y) >:= ∫ Y f (x, y)dy with Y the unit cell. Moreover, the effective density reads as One notes that the effective density is thus scalar valued and cannot achieve a rank-2 symmetric density tensor as required by the approximate Willis medium constituting the cloak, according to (7). Furthermore, it is suggested in 19 that at least in some cases, in order for one to perform cloaking with an elastic material defined by a symmetric elasticity tensor, one might need that the material be also defined by an anisotropic density. Inspired by earlier work on locally resonant acoustic metamaterials for anisotropic effective density, [28][29][30][31][32][33][34][35] we therefore add a second step in the homogenization process, that amounts to adding concrete bars in the soil. As is well known using spring-mass approximations of locally resonant elastic metamaterials, a dynamic effective (diagonal) density tensor is given by: 35 which is a dynamic correction to (14) with ∆ i the oscillator strengths of the resonances Ω i . Equipped with (13) and (15) for the effective anisotropic elasticity and density tensors, we can now design approximate seismic cloaks with concrete bars adding the required resonant feature to the soil.

V. ELASTIC THICK PLATES STRUCTURED WITH CONCRETE COLUMNS FOR HOMOGENIZED WILLIS MEDIUM WITHOUT RANK-3 TENSORS BUT ANISOTROPIC EFFECTIVE DYNAMIC DENSITY
In this section we investigate the case of seimsic cloaks which are a mitigation of cylindrical elastodynamic cloaks as proposed in Ref. 36 for in-plane coupled shear and pressure waves (which unfortunately require Cosserat media) and thin plate cloaks [40][41][42][43][44] for Lamb waves. Our case is indeed dedicated to Rayleigh like waves in thick plates. We do not claim that our cloaks work for body waves.
Bearing in mind that many authors have thus far achieved dynamic anisotropic density, see for instance, [30][31][32]34,35 symmetrized Willis cloaks are then achievable by classical homogenization approaches (getting rid of rank-3 tensors) with high contrast in material parameters for the artificial anisotropic density.
We would like to explore the propagation of elastodynamic waves in two types of cloaks that are in line with such a high-contrast homogenization, using the finite element method. The two cloaks are dedicated to seismic protection and could be achieved by manufacturing soil using concrete columns. Physical parameters are λ 0 = 87.7 * 10 6 Pa, µ 0 = 58.5 * 10 6 Pa and density ρ 0 = 1800 kg/m 3 for soil and λ 1 = 18.24 * 10 9 Pa, µ 1 = 9.39 * 10 9 Pa and density ρ 1 = 2300 kg/m 3 for concrete.
In order to prevent eventual reflections on the boundaries, Perfectly Matched Layers (PMLs) are applied to simulate the infinite extent of the elastic medium in the horizontal plane (note that we use socalled adaptative PMLs that work from high-frequencies down to the quasi-static limit, see Ref. 19). A point source is placed atop a semi infinite medium representing the earthquake epicenter, which provoke a Rayleigh-like wave (a type of wave that can be harmful and detrimental for buildings). The outside and inside radii of the cloak measure 21.6 m and 10 m, respectively. The point source is located at 50 m from the cloak center. Concrete columns are 2 m in diameter and penetrate 30 m in depth. Since we are interested in a frequency range greater than 5 Hz, the soil depth is set to 50 m (note that at 5 Hz Rayleigh waves have a wavelength of 36 m) and we apply free traction boundary condition to the bottom). We have checked that such a condition gives comparable results to those that are obtained if the bedrock is identified 50 m beneath soil. Indeed, potential reflection occur in both cases due to impedance mismatch. We deliberately consider such configuration as it seems more realistic than using PMLs that would model an infinitely deep, unstratified, soil. Concrete columns are arranged in a way that the effective density decreases as we move towards the center of the cloak. The designs of structured cloaks are optimized in order to achieve a reasonable balance between approximation of ideal cloak parameters according to effective equations of section IV and what can be implemented in practice with civil engineering technology at hand (i.e. various distributions of concrete columns of various diameter and length have been numerically tested). Indeed, the range of frequencies from 5 Hz to 10 Hz corresponds to wavelengths much larger than the columns's diameters (subwavelength inclusions underpins homogenization theory) and small enough compared to the basin depth to investigate Rayleigh waves rather than plate waves. Since the first cloak design, which we propose is deduced from elements of thin plate theory, we also look at the frequencies from 3 Hz to 5 Hz, which correspond to a transition between plate and Rayleigh waves. Indeed, the frequency range of interest for civil engineering is from 1 to 10 Hz, and it is interesting to check whether or not our cloaks's designs work beyond their 'legitimate' range of frequencies.

A. A first cloak design based on thin plate approximation
Nevertheless, we notice that the distribution of inclusions in Fig. 3(a)-(b) is coherent with the formula of the density of an elastic cloak previously reported by Ref. 40 in the limit of thin plates (so-called Kirchhoff-Love theory) and reads: with α = r 2 r 2 −r 1 . However, this step doesn't achieve the required effective anisotropy, and should be only valid for plate waves (which are akin to Rayleigh waves in the limit of thick plate not covered by Kirchhoff-Love theory).  (16) for cloak of type I (with a decreasing effective density towards the inner boundary of the cloak, as can be seen from the spatial distribution of inclusions in (a)-(b)) and of parameters in (6)-(7) for cloak of type II, according to effective equations (13)- (14), with an increasing anisotropy towards the inner boundary of the cloak (the inner boundary would approximate an infinite anisotropic tensor of elasticity in the limit of densely packed inclusions).
We depict in Fig. 4 the real part of the out of plane displacement field in 3D and top view representations at 5.5 Hz. It is observed that more than 50% of the field is suppressed at the center of the cloak (zone of protection) but recovered at its exit. The acceleration phenomenon is also noticed within the cloak area due to the higher effective Young's modulus. However, when we compute the cloak efficiency for protection throughout the frequency range [3Hz,10Hz], we FIG. 4. Mesh of computational domain (a) with 2,308,656 degrees of freedom, 562,675 elements and three-dimensional views (b,c) and top views (d,e) of the out of plane displacement field (Re(u z )) at 5.5 Hz of a Rayleigh-like wave propagating atop a soil without inclusions (b,d) and atop a soil manufactured with cloak of type I (c,e). One notes that the surface wave wavefront is nearly unperturbed by the cloak, but its amplitude is decreased within the center of the cloak. Fig. 4 with soil structured by a cloak of type II. Note the improved protection in the center of the cloak and the improved restitution of the field on the exit of the cloak compared with cloak of type I in Fig. 4. find that protection is merely 40% on average, see Fig. 7(b). This is quite a dramatic reduction of protection compared with the ideal cloak, which has 70% of protection on average, see Fig. 7(a).

FIG. 5. Same as
The cloak is efficient between 3 Hz and 10 Hz as it transpires in Fig. 7. It is worth noting that the lower limit (3 Hz) is connected to the basin depth (in which case waves transition from Rayeligh to Lamb) and the upper limit (10 Hz) to the diameter of the columns. The smaller the diameter the higher the upper frequency limit. Indeed the phenomenon is constrained by the diffraction that might be caused by the columns as can be observed in Fig. 6(a4).

B. A second cloak design based on approximate Willis equations in thick plates
We now make use of the much more involved transformational (section III) and effective (section IV) models for Willis equations in thick plates. In the configuration of cloak of type II, the inward effective anisotropy increases. This is accompanied by an undesirable reduction of the effective density the outer boundary of the cloak. This can be explained as the following: the concrete columns are arranged by bunch of clusters with one element over the first layer, two elements over the second and so forth up to the sixth layer. The seventh and last layer is isotropic and could be seen as a seismic shield. We depict in Fig. 5 the real part of the out of plane displacement field in 3D and top view representations at 5.5 Hz. For the sake of comparison, the free medium is also represented. We can notice that unlike the first configuration, there is much less reflection at the cloak entrance (due to the smooth mismatch of impedance), a quasi perfect restitution of the field at the exit and a FIG. 6. Top views of out of plane displacement field (Re(u z )) at 5.5 Hz (a1,b1), 7 Hz (a2,b2), 8 Hz (a3,b3) and 10 Hz (a4,b4). One notes that the wave wavefront is nearly unperturbed by the cloak of type II, and its amplitude is lower within the center of the cloak II than in the center of cloak I. comparable protection inside the cloak (almost 50% of the energy is suppressed on average in the frequency interval [5Hz,10Hz] as can be seen in Fig. 7(c)). However, it is noticed in Fig. 7(a), that some resonances occur in the ideal cloak, and similar effects occur with cloak of type II (notably around 8.5 Hz) in Fig. 7(c).

C. First conclusions on cloaks of type I and II
As can be seen in Fig. 6, cloak of type II performs clearly better than cloak of type I for both protection and invisibility. This can be attributed to the fact that the latter has been designed using transformed equations, 42 and homogenization techniques 41 in the context of thin plate theory. Kirchhoff-Love theory is an approach well suited for control of Lamb waves as experimentally validated in Ref. 43. However, such a limit theory is not legitimate in our case of thick plates, at least from a mathematical standpoint (the Rayleigh-like wave wavelength should be at least ten times smaller than the plate thickness to satisfy the hypothesis of thin plate, which is clearly not the case in our study). Interestingly, both cloaks of type I and II display resonances similarly to ideal cloak, see Fig. 7 (with a better match between ideal cloak and cloak of type II). We believe such unwanted resonances require a specific analysis of its own, as addressed in the next section. It is interesting to note that the level of protection displayed by our three cloaks (both ideal cloak and cloaks of types I and II) is lower than what is achieved through shielding effects via low frequency stop bands in Ref. 51. Although we consider the same frequency range as in Ref. 51, in our case we detour Rayleigh waves around the protected area, and there is little backscattered wave, so the protection mechanism is radically different.

VI. ESTIMATE OF LOCAL RESONANCES INSIDE THE SEISMIC CLOAK
Following, 26 where the case of Cosserat spherical cloaks is dealt with, in this section we investigate the origin of the resonances in the ideal cylindrical cloak by simplifying the governing Navier equations into a simple spring-mass model approximation leading to a transcendental equation. To do so, we assume that the plate thickness goes to infinity, in which case one can decouple the Navier equations as a vector problem for in-plane shear and pressure waves and a scalar out-of-plane problem for the vertical displacement u z . By doing so, we do not exactly address the problem of the Rayleigh waves, but this will nonetheless give us a reasonable estimate of the resonances, as we shall see. Because of the symmetry of the cloak, one can assume the out-of-plane displacement has a form u z =u z (r ), solution of the scalar Helmholtz equation where µ r (r ) = µ 0 (r − r 1 )/r , µ θ (r ) = µ 0 r /(r − r 1 ) and ρ z (r ) = ρ 0 ((r − r 1 )/r )(r 2 2 /(r 2 − r 1 ) 2 ). Actually, (17) can be written as a Bessel equation so it has a general solution, which is expressed in terms of cylindrical Bessel functions of first and second kinds (already derived in a different context in Ref. 45): Then, assuming that du z (r 1 )/dr=0 and u z (r 2 )=0, we find that the eigenfrequency should satisfy the following equation Y 0 (ωr 2 µ 0 /ρ 0 )/J 0 (ωr 2 µ 0 /ρ 0 ) = 0, which gives the frequency estimate f = ω/(2π) ∼ 6 Hz since the first zero of Y 0 (x) is 3.8317. This frequency estimate is in good agreement with the moderate elastic field enhancement shown in Fig. 7 around 5.3 Hz. The frequency estimate for the second resonance is f = ω/(2π) ∼ 11 Hz since the second zero of Y 0 (x) is 7.0156 and this is also in qualitatively good agreement with the second peak in Fig. 7 around 7.2 Hz. In fact, the predicted frequencies are overestimated since they are associated with out-of-plane shear waves, that propagate faster than Rayleigh waves. FIG. 7. Wave protection performances of the ideal cloaks (a) and approximate cloaks of types I (b) and II (c). Red (resp. green) curves correspond to integral of out-of-plane displacement computed over the surface (resp. over the volume) of the invisibility region (i.e. center of cloak) normalized by same in free space. Note that Rayleigh wave wavelength is 60 m at 3 Hz and 36 m at 5 Hz, whereas plate thickness is 50 m. Therefore, the frequency range [3Hz,5Hz] corresponds to a transition between plate and Rayleigh-like waves. Note also the strong similarities between resonances in (a) and (c); Indeed, cloak of type II in (c) is a better approximation of ideal cloak in (a), than cloak of type I in (b), the latter being designed using thin plate theory. Achieved protection is 70 on average in (a), 40 in (b) and 50 in (c).

VII. CONCLUDING REMARKS
In this article, we have conclusively shown that transformed Navier equations of the Willis type with a symmetric transformed elasticity tensor, allows for approximate cloak's designs getting rid of the rank-3 transformed tensors unveiled in Ref. 15. A further approximation consists in structuring soil with buried concrete pillars judiciously placed (note that this is the opposite case to what is proposed in Ref. 21, wherein soil was structured with softer columns). This is done and finite element computations confirm the cloaking efficiency (invisibility and protection) for two types of large scale seismic cloaks within a soft elastic plate (with soil parameters) 50 m deep reinforced with judiciously placed columns of concrete 20 m long and 2 m in diameter. Indeed, the two designs which we propose allow one to considerably reduce the elastic field vibrations in the center of the cloak with virtually no disturbance of the wave field outside the cloak. The protection could be further improved using the concept of mixed cloak, 26 which amounts to adding a PML layer at the inner boundary of the cloak, in essence some absorptive anisotropic medium. Of course, alternative routes exist to seismic wave protection, and we would like to mention very promising designs of seismic metamaterials in the tracks of Ménard's 2012 field test experiment in Grenoble, 14 notably so-called seismic waveguides, 46 but also large scale isochronous mechanical oscillators, 47 inertial resonators 4,48,49 and auxetic-like foundations 50 allowing for low frequency stop bands. Actually, in Ref. 51, numerical analysis of both surface and guided waves has been performed in miscellaneous large scale phononic crystals and metamaterials, with soil dissipation effects in viscoelastic media. Protection has been shown for frequency ranges associated with low frequency stop bands below 10 Hz. Such design of seismic shields could be combined with seismic cloaks to offer a wave protection for both surface and bulk waves which does not have adverse effects for the surrounding environment. Speaking of which the idea of forest of trees for seismic wave mitigation remains the most environmental friendly proposal, 52 although it has been observed thus far that wave protection is achieved beyond the 1 to 10 Hz frequency range of interest for civil engineers.
We would like to add that our results could be translated into geophysics upon use of time domain computations. It would be also important to study the effect of viscoelasticity on the elastic wave propagation (in a way similar to what was done in Ref. 51), which can be done with comsol multiphysics. Finally, some implementation of the seismic cloak with subsequent field test is in progress with the Ménard company. We show in Figure 8 a typical field test experiment by the Ménard company that could serve as a basis for the validation of our concept of seismic cloak, although this requires scaling down the diameter of inclusions of our study from 2 m to 0.4 m, as well as taking into account other constraints on minimum distance between inclusions, field topography etc. New designs of seismic cloaks are therefore underway to find a good compromise between cloak efficiency and feasibility in practice.