Liyun Cao

Zhichun Yang
email: yangzc@nwpu.edu.cn

Yanlong Xu

Shi-Wang Fan

Yifan Zhu

Zhaolin Chen

Yong Li

Badreddine Assouar
email: badreddine.assouar@univ-lorraine.fr

Flexural wave absorption by lossy gradient elastic metasurface

In this research, we systematically study the flexural waves diffraction. Based on the diffraction mechanism, we propose the concept of subwavelength lossy gradient elastic metasurface for flexural waves absorption. We theoretically reveal that the highefficiency absorption behavior stems from maximum multireflection-enhanced absorption of 0 th order diffraction, and experimentally show that robust flexural wave quasi-omnidirectional absorption in the frequency range extending approximately from 340 Hz to 1000 Hz (larger than 1.5 octaves). In addition, we propose a general approach which involves new physics of adjusting the arrangement sequence of subunits to suppress the 1 st diffraction mode, to further reduce the sub-wavelength thickness of the metasurface while maintaining its high-efficiency absorption. Our designs could provide new routes to broadband vibration suppression and cancellation in lowfrequency by lossy elastic metamaterials and metasurfaces.

Introduction

Vibration suppression of plate-like structures is of great significance in many engineering applications, for example, preventing vibration hazards and reducing noise.

Since early 1950s, this research topic has been widely concerned. The conventional suppression methods can be mainly divided into passive and active suppression approaches. The passive suppression method [START_REF] Kerwin | Damping of Flexural Waves by a Constrained Viscoelastic Layer[END_REF]Sun et al., 1995;Warburton, 1982), which is represented by attaching a large number of damping absorbers on the plate-like structures, is difficult to fulfil the requirements of lightweight and broadband features. The active suppression method [START_REF] Agnes | Active/passive piezoelectric vibration suppression[END_REF][START_REF] Dubay | Finite element based model predictive control for active vibration suppression of a one-link flexible manipulator[END_REF][START_REF] Moheimani | A survey of recent innovations in vibration damping and control using shunted piezoelectric transducers[END_REF][START_REF] Niederberger | An autonomous shunt circuit for vibration damping[END_REF], which is represented by applying control force through the actuator to the host structure based on feedback response, has complex structure and low stability. As a new generated method of vibration suppression, employing an acoustic black hole (ABH) [START_REF] Krylov | Experimental investigation of the acoustic black hole effect for flexural waves in tapered plates[END_REF][START_REF] Ma | Topological optimization of damping layout for minimized sound radiation of an acoustic black hole plate[END_REF][START_REF] Mccormick | Design optimization and performance comparison of three styles of one-dimensional acoustic black hole vibration absorbers[END_REF][START_REF] Pelat | Acoustic Focusing and Energy Confinement Based on Multilateral Metasurfaces[END_REF]Tang et al., 2016;Warburton, 1982), which can absorb the propagating waves (dominated by flexural waves) by little attached damping, promotes development of lightweight suppression technique. However, the ABH has weak damping effect at low frequencies and is not useful for a plethora of technological applications where the low frequencies are the main source to be damped. As another new generation of vibration suppression method, laying artificial periodic structures (Badreddine [START_REF] Badreddine Assouar | Broadband plate-type acoustic metamaterial for low-frequency sound attenuation[END_REF][START_REF] Fang | Ultra-low and ultra-broad-band nonlinear acoustic metamaterials[END_REF][START_REF] Oudich | A sonic band gap based on the locally resonant phononic plates with stubs[END_REF] on the plate-like structures can produce bandgap to prohibit wave propagation in low-frequency. Artificial periodic structures just isolate the vibration energy and do not dissipate the energy fundamentally. In addition, only subunit size is sub-wavelength, and the overall size of the periodic structure is still much larger than the wavelength. Recently, some vibration suppression techniques based on different methods, for example an open lossy resonators in onedimensional elastic beams [START_REF] Leng | Limits of flexural wave absorption by open lossy resonators: reflection and transmission problems[END_REF] and total reflection of propagating waves (Zhu et al., 2018a), have also been proposed. However, narrowband is their common characteristic. To date, it is still a great challenge to realize sub-wavelength broadband vibration suppression in low-frequency. Therefore, it is necessary to explore new mechanisms and approaches to break through it.

As one ultra-thin kind of metamaterials [START_REF] Dong | Robust 2D/3D multi-polar acoustic metamaterials with broadband double negativity[END_REF][START_REF] Kweun | Transmodal Fabry-Perot Resonance: Theory and Realization with Elastic Metamaterials[END_REF]Sugino et al., 2018;[START_REF] Zhang | An asymmetric elastic metamaterial model for elastic wave cloaking[END_REF]Zhu et al., 2014a, b), acoustic metasurfaces with sub-wavelength [START_REF] Assouar | Acoustic metasurfaces[END_REF][START_REF] Fan | Tunable Broadband Reflective Acoustic Metasurface[END_REF][START_REF] Li | Tunable asymmetric transmission via Lossy acoustic metasurfaces[END_REF]Qi et al., 2017;Zhu et al., 2018b) opened up new possibilities to realize extraordinary wave manipulation based on the generalized Snell's law (GSL) [START_REF] Yu | Light Propagation with Phase Discontinuities: Generalized Laws of Reflection and Refraction[END_REF]. However, with the development of research, it was found that the GSL will no longer be valid for some new phenomena involved high order diffractions, such as incident waves beyond the so-called critical angle (Liu et al., 2017a;[START_REF] Zhou | Three-dimensional acoustic characteristic study of porous metasurface[END_REF] and unexpected leakage waves (Cao et al., 2018b;Xie et al., 2014). For the high order diffractions, the directions of diffraction channels need to be predicted by the general formulas, i.e., the diffraction theorem [START_REF] Larouche | Reconciliation of generalized refraction with diffraction theory[END_REF]. Normally, several diffraction channels simultaneously exist for a particular incidence, and the magnitudes of diffractions in different channels will be different. Apparently, the complicated diffraction will not be fully revealed by only the directions of the diffraction channels. The magnitude of each order diffraction should also be explored. For revealing the comprehensive underlying diffraction mechanism, the mode-coupling method [START_REF] Liu | Metasurface-based angle-selective multichannel acoustic refractor[END_REF]Mei and Wu, 2014;[START_REF] Hou | Highly Efficient Acoustic Metagrating with Strongly Coupled Surface Grooves[END_REF] had been used to accurately solve the amplitude of each order diffraction. Recently, based on these methods, the acoustic diffraction mechanism had been studied systematically and some new physics related to acoustic diffraction had been revealed. As typical representatives, acoustic asymmetric transmission [START_REF] Li | Tunable asymmetric transmission via Lossy acoustic metasurfaces[END_REF] and acoustic absorption [START_REF] Shen | Harnessing Multiple Internal Reflections to Design Highly Absorptive Acoustic Metasurfaces[END_REF]Yi Fang, 2018;[START_REF] Zhou | Three-dimensional acoustic characteristic study of porous metasurface[END_REF] can be obtained by high order diffraction. Anomalous transmission and reflection through high order diffraction can be completely reversed by changing the integer parity of the phase gradient metagratings [START_REF] Fu | Reversal of transmission and reflection based on acoustic metagratings with integer parity design[END_REF].

Recently, although big progresses have also been subsequently achieved in the field of elastic metasurfaces (Liu et al., 2017b;Xu et al., 2019b;Cao et al., 2020a;Cao et al., 2018a;[START_REF] Chen | A programmable metasurface for real time control of broadband elastic rays[END_REF][START_REF] Lee | Mass-stiffness substructuring of an elastic metasurface for full transmission beam steering[END_REF][START_REF] Li | Tunable modulation of refracted lamb wave front facilitated by adaptive elastic metasurfaces[END_REF]Liu et al., 2017b;[START_REF] Park | Recent Advances in Non-Traditional Elastic Wave Manipulation by Macroscopic Artificial Structures[END_REF][START_REF] Qiu | Steering and focusing of fundamental shear horizontal guided waves in plates by using multiple-strip metasurfaces[END_REF][START_REF] Zhang | Metasurface constituted by thin composite beams to steer flexural waves in thin plates[END_REF]Tian and Yu, 2019;Xia et al., 2019;Xu et al., 2019a;[START_REF] Yuan | Switchable multifunctional fish-bone elastic metasurface for transmitted plate wave modulation[END_REF][START_REF] Zhang | Metasurface constituted by thin composite beams to steer flexural waves in thin plates[END_REF][START_REF] Zhu | Anomalous Refraction of Acoustic Guided Waves in Solids with Geometrically Tapered Metasurfaces[END_REF], almost all existing researches on elastic metasurfaces are simply related to the manipulation of the wavefield based on the GSL, which is only related to 1 st order diffraction. The manipulation of other orders of elastic wave diffractions based on elastic metasurface has not been systematically studied due to the essential difference between the elastic waves and the acoustic waves; for example, more degrees of freedom and more complexity of the control equation. We expect to develop an effective theoretical tool, which is similar to the mode-coupling method in acoustics (Mei and Wu, 2014;[START_REF] Hou | Highly Efficient Acoustic Metagrating with Strongly Coupled Surface Grooves[END_REF], and explore the rich physics of higher-order diffractions in the elastic metasurface. Fortunately, similar methods (Willis, 2016;Srivastava, 2016), which are also called mode-coupling methods for consistency, had been used in elastic waves to study negative refraction of anti-plane shear waves at a plane interface between a homogeneous elastic half-space and a layered periodic composite. Recently, these methods [START_REF] Lustig | Anomalous energy transport in laminates with exceptional points[END_REF][START_REF] Lustig | On the band gap universality of multiphase laminates and its applications[END_REF]; [START_REF] Mokhtari | On the properties of phononic eigenvalue problems[END_REF][START_REF] Mokhtari | Scattering of in-plane elastic waves at metamaterial interfaces[END_REF] had been extended to study the scattering of in-plane elastic waves.

Although they only focus on anti-plane shear waves or in-plane waves, they provide a good theoretical basis for the studies of other elastic wave patterns in different periodic structures.

In the present study, we expand the conventional mode-coupling method to systematically study the diffraction mechanism of flexural waves. Based on the diffraction mechanism, we introduce the concept of lossy gradient elastic metasurface (LGEM) to explore the underlying physics of flexural wave absorption. Different from the inherent energy loss in acoustic waves, the small damping property of the solid material makes the loss of elastic wave in the structure negligible. Therefore, we introduce different lossy physical system into the elastic metasurface, which is an additional constrained damping layer (i.e., a damping layer of butyl rubber and a constraint layer of aluminum foil) on the subunits. For the lossy subunits, the equivalent models are established to analytically predict the amplitudes and phases of reflected waves. Further, we theoretically analyze the absorption performance in details. In addition, we propose a general approach to reduce the sub-wavelength thickness of the metasurface while maintaining its efficient absorption. Finally, both numerical simulations and experiments are carried out to demonstrate the broadband and highefficiency flexural waves absorption of the LGEM in low-frequency.

Design of the subunits of the LGEM

Effective model of the lossy unit

First, the subunits without the damping and constrained layers are studied, i.e., d1 = d2 = 0. The one-dimensional governing equation for the flexural wave in the plates of Regions (I) and (II) can be expressed as the following form: ,

where is Young's modulus of the plate, is the density of the plate and is the moment of area, in which and are the thickness and

Poisson's ratio of the plate, respectively.

The one-dimensional governing equation of flexural waves is a fourth-order partial derivative equation. The wavenumber has four solutions, i.e., two real wavenumbers ( ) ( ) Therefore, the general solution of displacement for the governing equation is:

4 2 0 0 0 0 4 2 , , 0 w x t w x t E I d x t r ¶ ¶ + = ¶ ¶ 0 E 0 r 3
,

where , , , and are complex coefficients. and correspond to the positive-going and negative-going propagating flexural waves, respectively. and correspond to the positive-going and negativegoing evanescent flexural waves, respectively. The real wavenumber can be expressed as , in which is the circular frequency.

Further, the subunits with damping and constrained layers are studied, as shown in Fig. 1(b). To simplify, the lossy subunit with a three-layer composite structure can be considered as an effective one, i.e., an isotropic plate, as shown in Fig. 1(c). The effective bending stiffness of the lossy unit can be expressed as [START_REF] Ross | Damping of plate flexural vibration by means of viscoelastic laminate[END_REF]:

.

(3)

Based on the effective bending stiffness, we can obtain the effective wave number in the lossy unit, which is corresponding to the wavenumber of the positive-going propagating mode in Eq. ( 1). The detailed solutions of the effective bending stiffness and effective wave number are reported in Appendix A. Due to the loss of the damping layer, the solved effective wave number has a small imaginary part. It can be expressed as , where and are the real and imaginary parts,

0 k ± 0 i k ± × ( ) (
)

0 0 0 0 - - 0 0 0 0 0 , ik x ik x k x k x i t w x t A e B e C e D e e w = + + + 0 A 0 B 0 C 0 D 0 0 ik x A e - 0 0 ik x B e 0 0 k x C e - 0 0 k x D e 0 k 1/4 2 0 0 0 0 0 d k E I r w ae ö = ç ÷ è ø 2 f w p = ( ) ( ) ( ) ( ) 2 2 2 eff 0 1 2 0 1 20 2 30 20 30 1 1 23 0 + + + 12 2 1 c c c c c c c c c c D D D D K d K d d K d d d d d d d K K d d g é ù = + - + - ë û é - ù - ae ö - + + - ç ÷ ê ú + è ø ë û eff k Re Im eff eff eff k k i k = -× Re eff k Im eff
k respectively. Furthermore, the displacement of positive-going propagating flexural wave in the lossy subunit can be rewritten as , where and represent fluctuation harmonically and decay exponentially in space, respectively. Therefore, and correspond to the amplitude and phase shift of the reflection wave emanating from the lossy subunit, respectively.

In order to obtain the effective mechanical properties of the lossy subunit, we measure the storage modulus (the real part of the complex Young modulus) and the loss factor ƞ of its damping layer varying with the frequency by the Dynamical Mechanical Analysis (DMA). The test set-up is shown in Fig. 2(a), and the test sample of butyl rubber is shown in the illustration in the lower left corner. The operating frequency range of our DMA is from 0 to 1000 Hz. The following investigations will be based on the measured parameters in this frequency range, which is sufficient for the following theoretical analysis of vibration absorption. The measured datum of storage modulus and loss factor are shown in Fig. 2(b). For convenience, the fitting curves of the datum are obtained, and the corresponding fitting functions can be expressed as . ( 4)

The influence of constrained and damping layers parameters on the subunit

Without considering the constrained and damping layers, the phase shifts of the reflected waves emanating from the j th subunit can be simply

E f f f f f f f f f h - - - - = × × - × × + × × £ £ - × × + × × - × × + × 3 10 0.225, ( 200
Hz 1000 Hz) f f - ì ï ï ï í ï ï ï × + £ £ î 0 2 j j k h f =
controlled by the lengths . When these lengths of the subunits fulfil (Cao et al., 2018b) ,

the phase shift of the reflected wave on the surface of the gradient elastic metasurface without a loss (GEM) is approximately linear, where and are the wavelength and additional fixed length, respectively. In this way, the reflected wave can be manipulated by the GEM based on the GSL and diffraction theorem. In these subunits, the phase shift between the two subunits of the 1 st subunit and J th subunit is the maximum, i.e., . After introducing the constraint damping layer to the subunits, the maximum phase shift will be rewritten as .

For independently defining the effect of loss and phase shift on absorption, when the difference between the phase resolution for lossless and lossy subunits is less than about rad, the phase difference caused by the constraint damping layer can be ignored. Therefore, the difference between the two maximum phase shifts needs to fulfil .

In this way, the phase gradient designed by Eq. ( 5) can still be approximately linear for the lossy subunits. According to Eqs. ( 5)-( 6), we can get .

j h 0 0 , ( =1,2 ) 2 j j h h j J J l = + ! 0 0 =2 / k l p 0 h 1 0 Max( )=2( -) J h h k f Re loss 1 eff Max( ) =2( -) J h h k f 0.32% 2 0.02 p × » ( ) loss = Max( )-Max( ) 0.02 -1 rad J f f f D < FIG. 2 (a)
(7)

The central frequency in our design is 600 Hz, corresponding to a wavelength of . Without loss of generality, the number of subunits and additional fixed length are chosen as J =12 and , respectively. In order to obtain a quantitative evaluation of the dependence of the real and imaginary parts of on the thicknesses of the damping layer d1 and the constrained layer d2, we define the correlation coefficient as Recalling Eq. ( 7), we can get where the phase difference caused by the constraint damping layer can be ignored. Therefore, the corresponding the correlation can be chosen as 0.9. According to Figs. 2(c) and 2(d), the thicknesses of the damping layer and the constrained layer is 3 mm and 0.15 mm, respectively.

Re eff k D Im eff k eff k Re Re eff eff 0 0 0.02 = J k k k l D -< 0 =221.2 mm l 0 0 / 2 h l = eff k , ( 8

The amplitude and phase shift of the subunits

For the subunits with the certain thicknesses of constrained and damping layers, and varying with frequency are shown in Fig. 3(a), respectively. It can be seen that the value of is less than 1.08 in the whole frequency range and increases with the frequency. This verifies that phase difference caused by the designed constraint damping layer can be ignored. Further, we accurately solve the amplitude and phase shift of the reflected waves emanating from the lossless and lossy subunits. where is the position coordinate of the j th subunit. We further obtain the transfer equation between the state vectors of the right end in Region (I) and the left end in Region (II) by integrating Eq. ( 11) with dy at the region

f f k d d k d d f df f f D = D - ò ( ) ( ) max min Im Im eff 1 2 eff 1 2 max min 1 , , , f f k d d k d d f df f f = - ò min f max f Re eff k D Im eff k Re eff 0 0.02 / 1.08 k J l D < » Re eff k D Re eff k D Im eff k Re eff k D Im eff k w j V M A B C D { } T , , , w M V j = v { } , , ,
, ( 12
)
where l is the total width of subunit and slit, . According to Eqs. ( 10) and ( 12), we can get

j w j ¢ = = M EIw¢¢ V EIw¢¢¢ = - ¢ (I) (I) 1 (II) (II) 2 = = T T × × k k v v 1 T 2 T (I) (II) R L (I) (II) R L (II) (I) L R (II) (I) L R , , 0, other , 0, other j j w w M y y p M V y y p V j j = = ì < + = í î ì < + = í î , j y / 2 j y y l < + (II) (I) (I) L R 3 R 1 0 0 0 0 1 0 0 = 0 0 0 0 0 0 T e e é ù ê ú ê ú = ê ú ê ú ë û v v v l p e = . ( 13
)
It should be pointed out that when the transfer matrix of coefficients vectors is an identity matrix, the impedances of Regions (I) and (II) at the interface are matched. It means that the propagating wave will not be reflected from the interface.

On the other hand, it needs to fulfil the requirement of free boundary conditions in the right boundary of Region (II), i.e., the shear force and the bending moment should be zero. It gives

, ( 14
)
where is the matrix for the free boundary conditions and given in Appendix B. The wave field for Region (I) can be described as , where

is the incident flexural wave with an amplitude of 1. and are the amplitude ratios of the reflected propagating flexural wave and reflected evanescent flexural wave in Region (I), respectively. Thus, we can get the coefficient vectors . According to Eqs. ( 13) and ( 14), we have

, ( 15
)
where T5 is the transfer matrix for the j th subunit. The amplitude ratio and corresponding phase of the reflected propagating flexural wave in the far field of Region (I) can be solved by Eq. ( 15) (Cao et al., 2020b;Xu et al., 2019a). The phase shifts of the reflected waves for all subunits are obtained by subtracting their phase from the phase for the 1 st unit. When keff and Deff are equal to k0 and D0, respectively, the transfer matrix T5 of the lossy subunits will degenerate into that of the lossless subunits.

( ) ( ) The phase shifts of the lossy subunits are almost the same as that of the lossless subunits, while the amplitudes of the lossy subunits are reduced by about 80%. To quantify the performance of the subunits, we further examine the phase shifts and amplitudes

I I (II) 1 L 2 3 1 R t R T TT T - = = k k k t T (II) 4 L =0 T k 4 T ( ) 1 1 1 (I) * b b b ik x ik x k x b b w x e re r e - = + + 1 b ik x e - b r b r * (I) R = 1, ,0, T b b r r * é ù ë û k ( ) ( ) I I 4 R 5 R 0 t T T T = = k k FIG. 3 (a)

Re eff k D

Im eff k varying with the frequency for the lossy subunits (see Figs. 3(c) and3(e)) and lossless subunits (see Figs. 3(d) and3(f)). The phase shifts of the lossy subunits are almost the same as that of the lossless in the whole frequency range, while the amplitude of the lossy subunit decreases with frequency. This again verifies the phase shift caused by the loss in the designed subunits can be ignored. Furthermore, we can independently explore the effect of loss and phase shift on flexural wave absorption.

Analysis of flexural wave diffraction

Mode-coupling method for flexural waves

The above designed subunits are arranged periodically to compose the LGEM, as shown in Fig. 1(a). The reflection angles of these diffraction modes can be calculated by the diffraction theorem

, ( 16
)
where and mean y-component wave vectors of the incident waves and n th order diffraction, respectively. and are the reflection angle of the n th order diffraction and the incident angle. describes the phase gradient along the LGEM, and is the width of the supercell. When n = 1, Eq. ( 16) is transformed to , yielding the well-known GSL. According to Eq. ( 16), we can get the reflection angle of the n th order diffraction by . ( 17)

Further, we expand the mode-coupling method to calculate the reflection amplitudes of these flexural wave diffraction modes. In Fig. 1(a), the host plate and the LGEM

in yn y k k ng = + in 0 i n = sin y k k q 0 r = sin n yn k k q r n q in q 2 = g p g g in yn y k k g = + r in 0 sin sin n n k g q q = +
structure are divided into Region (P) and Region (S). For the sake of universality to study the vertically and obliquely incident flexural waves in the host plate, the governing equation for the flexural wave in Region (P) should be reconsidered in the two-dimension form:

, ( 18
)
where .

The displacement field including all diffraction modes in Region (P) can be expressed as:

, ( 19
)
where is the Kronecker delta, Ai is the amplitude of the incident wave, An and Bn are the reflection amplitudes of the n th order propagating and n th order evanescent flexural wave diffraction modes, respectively. and are x-component wave vectors of the propagation and evanescent flexural wave diffraction modes, respectively.

Recalling Eq. ( 17), when the phase gradient fulfils the relation of , we always can find an incident angle in the range from -90° to 90°

to make the reflection angle of the th order diffraction exist, where is the maximum of all the existing orders. In other words, there will exist the propagating diffraction modes of from the th order to the th order when the phase gradient fulfils the relation of , while other order diffraction modes are evanescent.

( ) ( )

2 P 2 2 0 00 2 + , , 0 D d w xyt t r ae ö ¶ Ñ Ñ = ç ÷ ¶ è ø 2 2 2 2 2 = x y ¶ ¶ Ñ + ¶ ¶ ( ) 0 0 P ,0 0 ( , ) ( ) y y n y n x x n x n ik y ik y ik y ik x ik x k x n i n n n w x y Ae e A e e B e e d ±¥ - × - × - × - × × × = = × + + å ,0 n d 2 2 0 = xn yn k k k - 2 2 0 ˆ= xn yn k i k k --- g 0 2 n N n k g = < in q rN q N N N - N g 0 2 ˆk N g <
The infinite summation in Eq. ( 19) includes only infinite propagating diffraction modes. The number of propagation modes can be determined by , (20) where "roundup" is a function which rounds up to the nearest integer. Therefore, the total number of those propagating modes is controlled by the phase gradient . The coefficient vector of the reflection diffraction field in Region (P) can be defined as , where .

In Region (S), since the width of the subunit is much smaller than the operating wavelength, only the fundamental mode needs to be considered. The displacement of flexural waves in the j th subunit can be expressed as ,

where , , , and are the amplitude coefficients. The coefficient vector of the fundamental mode in Region (S) can be defined as , where .

In the direction perpendicular to the periodic waveguide plates, waveshapes related to the y-component wave vectors of the reflected diffraction modes have the orthogonal relation:

, ( 22
)
where . Continuity (11) and orthogonality ( 22) are employed and modecoupling method is expanded to calculate these flexural wave diffraction modes. We get a linear equation set about and :

( )

0 =2 +1=2 roundup 2 / 1 N N k g × - g [ ] P r 1 = , , T n N Τ e e e ! ! [ ] , n n n A B = e eff eff eff eff S ( ) ik x k x ik x k x j j j j j w x a e b e c e d e - × -× × × = + + + ( ) j a j b j c j d 1 = , , T S j J é ù ë û Τ k k k ! ! = , , , j j j j j a b c d é ù ë û k ( ) 2 2 0, g j n g dy j n y y * - × = ¹ ò yn ik y n e y - × = P r Τ S Τ , (23)
where is a column vector of size whose elements are all Ai, and are square matrices of size . The detailed derivation can be found in Appendix C. The corresponding amplitudes of diffraction modes can be calculated by solving Eq. ( 23). The reflection coefficient of the n th order diffraction mode can be defined as

, ( 24
)
which represents the ratio of reflected energy to incident energy for n th order mode in

x-axis direction. When keff and Deff are equal to k0 and D0, respectively, Eqs. ( 23) and

(24) will degenerate to calculate reflection coefficient of diffraction modes for the lossless GEM.

The reflection angle and amplitude of each diffraction wave

First, without loss of generality, we begin with a lossless GEM with the phase gradient of . According to Eq. ( 20), three propagation modes can be predicted in the reflection field. The corresponding reflection angles of these modes are calculated by Eq. ( 17), as shown in Fig. 4(a-i). And the reflection coefficients associated with each order diffraction are obtained by Eqs. ( 23) and ( 24), as shown in Fig. 4(a-ii). The reflection coefficient associated with the 1 st or -1 th order mode is 1 at a relatively small incident angle, while the 0 th order mode dominates at a big incident angle. After the loss is considered in the LGEM, the corresponding reflection angles and reflection coefficients of these diffraction modes are shown in Figs. 4(b-i) and 4(b-ii). It can be 4(a-iii) and 4(b-iii) that the simulation results are in good agreement with the analytical ones. This confirms the accuracy of the proposed analytical method.

P P r i S ae ö × = × ç ÷ è ø 1 2 Τ G G Τ Τ P i Τ (4 2 ) J N + 1 G 2 G (4 2 ) (4 2 ) J N J N + ´+ 2 0 2 0 n n n x i k A r k A × = ×
It should be pointed out that the sum of the reflection coefficients of all different order modes for the lossless GEM is always 1 at different incident angles, as shown in Fig. 4 (a-ii). The analytical and simulated reflection wave fields in Fig. 4 (a-iii) also

show that the reflection amplitudes are 1. In other words, the incident energy of flexural waves is proved to be equal to the reflected energy. Therefore, there is no mode conversion in the GEM structure. In addition, since it is known that the loss cannot induce mode conversion, the mode conversion does not exist in the LGEM.

Design of the quasi-omnidirectional LGEM and the internal mechanism

Propagation diffraction modes induced by multiple reflections

In order to reveal the underlying physics of the disappearance of -1 th order diffraction mode after considering damping in Fig. 4(b-ii), we reexamine the diffraction theorem of Eq. ( 16). The physical meaning of in Eq. ( 16) is an extra momentum in y direction, where the value of n corresponds to n th order diffraction.

When n = 1, the extra momentum corresponds to 1 st order diffraction and the diffraction

0 =k g = 2 / n n g g p ×
theorem of Eq. ( 16) degenerates into the GSL. The total phase shift of the supercell consisting of J subunits is 2π, so the phase shifts between adjacent subunits is .

It is known for the GSL that the phase shifts of these subunits are caused by one-time reflection of the incident wave in the subunit. When , the total phase shift of the supercell is , which corresponds to the th order diffraction. The phase shifts between adjacent subunits increase to times, i.e., , while the geometric size of subunits remains unchanged. The phase shifts of these subunits need to be matched by -times reflection of the incident wave in the subunit. Therefore, th order diffraction corresponds to -times reflection in the subunit.

Similarly, when , the total phase shift of the supercell consisting of J subunits is . Therefore, the phase shifts between adjacent subunits is . But the phase shift is a negative value, which doesn't fit the actual physics.

Due to the periodicity of waves, the negative phase shift is equivalent to

. It needs to be matched by E-times reflection of the incident wave in the subunit. Recalling the phase shift between the adjacent subunits is for the one-time reflection of the incident wave in the subunit. Therefore, we can get

. ( 26
)
After simplification, Eq. ( 26) can be expressed as . Therefore, the time of multiple reflections in the subunit corresponding to the th order diffraction can be expressed as

2 J p 0 n N * = > 2 N p * × N * N * 2 N J p * × 2 N J p * × N * N * N * 0 n N * = £ 2 N p * × 2 N J p * × 2 + 2 N J p p * × 2 J p 2 2 2 + = N E J J p p p * × × E J N * = + N * . ( 27
)
According to Eq. ( 27), we can calculate the times of multiply reflection in the subunit for all propagation diffraction modes in Fig. 4(a-ii). The times of multiply reflection for the -1 th order, 0 th order, and 1 st order diffraction modes are 11, 12, and 1, respectively. On the other hand, the lower the corresponding times of multiple reflections for the existing propagation diffraction modes is, the less the resistance from the structure is. The diffraction mode with the least times of multiply reflections will dominate in all existing propagation modes. When the incident angle is the positive, existing propagation modes are the -1 th order and 0 th order ones. The times of multiply reflection for the -1 th order one is less than that for the 0 th order one. Therefore, the -1 th order diffraction dominates, as shown in Fig. 4(a-ii). Similarly, when the incident angle is negative, the 1 st order diffraction dominates. In Fig. 4(b-ii), when the damping is taken into account, 11-times multiply reflection for the -1 th order diffraction will enhance damping dissipation of flexural waves. It leads to that the absorption of the -1 th order diffraction is greater than that of the 1 st order diffraction. This is the reason for the disappearance of the -1 th order diffraction mode after considering damping.

Verification of multiple reflections

In order to verify the above multiple reflections, the case that the incident wave with the incident angle of 45 o propagates into the LGEM with the phase gradient of is adopted. The central frequency is 600 Hz. The model of the LGEM is shown in Fig. 5(a). The corresponding lengths of the subunits are , , ,

, 0 = , 0 N N E J N N * * * * ì > ï í + £ ï î 0 =k g 1 h 2 h ! -1 J h J h
respectively, which can be solved by Eq. ( 5). According to the above analysis, the -1 th order diffraction dominates in the reflection wave filed, which will be absorbed by the

LGEM. The absorption coefficients considering all diffraction modes can be defined as . We calculate the absorption coefficients varying with the number of subunits J in the supercell according to Eqs. ( 23) and ( 24). The results are shown in Fig. 5(c).

For comparison, we design another model, as shown in Fig. 5(b), in which all subunit lengths are set to the same value of . The model is named as the lossy uniform elastic metasurface (LUEM). The corresponding widths of the subunits and supercell are the same with the LGEM. When the uniform subunit length of LUEM fulfils , the damping area of the LUEM is the same as that of the LGEM.

If the incident wave is reflected by one time in the subunit of the LGEM, the absorption coefficients for the LUEM and the LGEM will be the same. However, in the section 4.1, it shows that the incident wave with the incident angle of 45 o will be reflected by J-1 times in the LGEM due to the dominant -1 th order diffraction.

Therefore, we increase the subunit length of the LUEM to . In this way, we analytically calculate the absorption coefficient of the LUEM for the incident wave with the incident angle of 45 o according to Eqs. ( 23) and ( 24). The results are depicted in Fig. 5(c). It can be seen that the absorption coefficient of the LUEM is exactly the same as that of the LGEM with different numbers of subunits. This result verifies that the times of multiple reflections for the -1 th order diffraction is J-1 calculated by Eq. LGEM and the LUEM.

n 0 =1 n r a ±¥ = -å h 1 1 J j j h h J = = å ( ) -1 J h × (27).

Quasi-omnidirectional flexural wave absorption of the LGEM

According to Eq. ( 27), we know that the number of multiple reflections for the 0 th order propagation mode in the LGEM is J times, which is the highest in that for all propagation modes. According to Eq. ( 20), we can adjust the surface phase gradient of the LGEM (J=12) as , so that there is only one propagation mode, i.e., the 0 th order one, in the reflection field. In the later research, we all adopt the LGEM with the phase gradient of . When the loss is not considered, the corresponding reflection angles and reflection coefficients of these diffraction modes are calculated analytically, as shown in Figs. 6(a-i) and 6(a-ii), respectively. It can be seen that the amplitudes of the 0 th order modes are always 1 for all incident angles in Fig. 6(a-ii).

When the loss is considered, the reflection coefficients jump from 1 to less than 0.2 for the incident angles from -75 o to 75 o due to the multiple reflections, as shown in ( )

1 J h -× 0 =2.1k g 0 =2.1k g Fig. 6(b-ii).
The LGEM realizes a high-efficiency and widen-angle absorption.

Therefore, it is a perfect candidate for quasi-omnidirectional flexural wave absorption.

To intuitively display the high-efficiency absorption, the incident wave with the incident angle of 0 o is chosen to propagate into the GEM and LGEM, respectively. For these two cases, full reflected wave fields are solved analytically and numerically, as Further, for the quasi-omnidirectional LGEM, we intuitively show the multiple 0 =2.1k g reflections in the subunits by the analyses of the time domain signal. The computational model without considering loss is shown in Fig. 7(a). For comparison, the model of the lossless UEM is also shown in Fig. 7(c). The only difference between the two models is the subunit length, the one is gradient and the other is uniform. The periodic boundary conditions are applied to the two long boundaries of the models.

Performance evaluation of the quasi-omnidirectional LGEM

To quantify the absorption performance of the quasi-omnidirectional LGEM, we examine the absorption coefficient varying with the frequency and incident angle by Eqs. ( 23) and ( 24). The corresponding results are shown in Fig. 8(a) with the white contour

. The area of the high absorption coefficient, which is defined when the absorption coefficient exceeds 0.8, is located in the frequency range extending approximately from 340 Hz to 1000 Hz and for angles between -75° to 75°, except for the small area in the upper left corner in Fig. 8(a). Next, we will reveal the reasons for the low-efficiency absorption region in Fig. 8(a), including the region in the upper left corner, large incident angle region and low-frequency region.

The low-efficiency absorption induced by 1 st order diffraction

In order to find out the reason of weak absorption in the upper left corner in Fig. 8(a), we calculate the reflection coefficients of the 0 th order, -1 th order, and 1 st order diffraction varying with the frequency and incident angle, without considering the damping. The results are shown in Figs. 8(d), 8(e), and 8(f), respectively. It can be seen that the -1 th order and 1 st order diffractions dominate in the upper right and left corners in these figures, respectively, while the 0 th order diffraction dominates in the rest of the region.

The reason is that when the frequency is greater than 660 Hz, the surface phase gradient will be less than . According to Eq. ( 20), the number of propagation modes increase to 3 for the incident angle range from -90° to 90°. In other words, another order propagation mode will exist in addition to the 0 th order diffraction. The added propagation modes are the -1 th order and 1 st order for the positive and negative incident angle, respectively, according to Eq. ( 17). On the other hand, we know from the above analysis that the 0 th order and -1 th order diffractions will be greatly absorbed due to multiple reflections while the 1 st order diffraction cannot do due to only one reflection. Therefore, in the upper left corner region of Fig. 8(a), the 1 st order diffraction dominates and its low-efficiency absorption makes the value of absorption coefficient very small. , where is the scattering matrix for a reciprocal system, which already contains the information of the evanescent wave. On the other hand, a general form of the 2 × 2 scattering matrix for a power-conserving and reciprocal system can be expressed as (Liu et al., 2017b)

Impedance mismatch at the interface

1 l { } 1 1 1 1 ˆˆ, , , T A B C D = 1 e { } 2 2 2 2 ˆˆ, , , T A B C D = 2 e 1 l 1 ˆˆ T = 1 2 e e { } { } 2 1 2 1 2 ˆˆˆ, = , T T A C T A C 2 T , ( 28
)
where has the physical meaning of power conservation. By comparing the two form of , , we can get the analytical impedance at the interface

. ( 29
)
It can be seen from the above expression that the impedance at the interface is independent of the frequency. Based on Eq. ( 29), the reflection coefficient can be obtained from and expressed as . interface . This indicates that for the incident wave with a big incident angle (>75°), the unabsorbed reflected wave from the LGEM is caused by the impedance mismatch at the interface between the host plate and the waveguide subunit. This is the reason that the absorptions coefficient sharply drops for the incident waves with very large incident angles. The phase shifts of all subunits for the frequencies of 343 Hz, 600 Hz, and 1000 Hz are shown in the illustration. (d), (e), and (f) are the reflection coefficients of 0 th order, -1 th order, and 1 st order diffractions varying with the frequency and incident angle, respectively, without considering the damping.

2 2 1 1 1 ˆ2 / 1 1 1 Z Z Z Z T Z Z Z Z h h ae ö - ç ÷ + + ç ÷ = ç ÷ - ç ÷ è + + ø h 2 T p l » 3 2 2 2 2 2 0 0 0 0 0 2 0 0 0 0 ( 1) ( 1 2( 1) 4 4 -) = y y x x y x k k k k k Z k k k k u u u u u - + + - - 1 ˆ1 Z r Z - = + 2 2 0 0 0 0 0 0 0 2 2 2 0 0 0 0 0 0 0 2 0 2 2 0 2( 1) 4 + ˆ( ) ( 1) ( ) ( 2( ) ) ) 1 1 ( 4 y x y x x y x y x x k k k k k k k r k k k k k k k k k u u u u u u - + + = - + + - - + - ) - ( 1

Degenerated phase gradient in the low-frequency

It is worth noting that in Fig. 8 The value of can quantify the degradation of the linear phase gradient. When is equal to 0, it shows that there is a perfect linear gradient. The more the value of is, the greater the degradation of the linear gradient is. According to Eqs. ( 15) and ( 31), the value of varying with the frequency is shown in Fig. 8(c). It can be seen that the change rate of has been less than -0.27 when the frequency is less than 343 Hz.

The phase shifts of all subunits for the frequencies of 343 Hz, 600 Hz, and 1000

Hz are shown in the illustrations of Fig. 8(c) as three typical examples. It can be seen that phase shifts of all subunits decrease in the frequency of 343 Hz, in other words, the linear phase gradient from the subunits has degenerated. It indicates that in the frequency range below 343 Hz, the 0 th order diffraction induced from the phase gradient has been greatly generated by the decreased phase shift of the subunits, which will not cause J-times reflections in the subunits. Therefore, the absorption coefficient sharply drops in the range below 343 Hz. It should be pointed out that in Fig. 3(c) the loss of the subunits has also some degeneration in the frequency range from 200 Hz to 343 Hz, but the reduced loss has little effect on the absorption coefficient. In addition, the degenerated phase gradient in the high frequency does not affect the 0 th order diffraction, which can result in high-efficiency absorption. The detailed explanation is shown in

Appendix D. ( ) ( ) 0 0 ˆ= f f f f d f - 0 ˆ=2 f p d d d d d

Reducing the thickness of the quasi-omnidirectional LGEM

The above LGEM is designed in the central frequency of 600 Hz and its thickness (in the propagation direction) is determined by the maximum subunit length of according to Eq. ( 5). When the frequency is less than 600 Hz, the designed thickness will be less than the wavelength. In other words, in the frequency range from 343 Hz to 600 Hz, high-efficiency absorption can be obtained by the LGEM with the sub-wavelength thickness, which is 0.76 times the maximum wavelength. In order to make the LGEM structure more compact in lower frequencies, we propose a general approach to further reduce its thickness while maintaining its high-efficiency absorption.

To understand, we first consider the supercell composed of 13 subunits, which is the 1 st case marked as "C1", shown in Fig. 9(a). The supercell has the same subunit width p and slit width as these of the above LGEM. All subunit lengths are redesigned to keep linear phase gradient by Eq. ( 5). According to Eq. ( 27), the time of multiple reflections for the 0 th order diffraction is increased to 13 due to the subunit number is changed from 12 to 13. The corresponding absorption coefficient is almost the same as that in Fig. 8(a) due to the slightly increased reflection time. Further, we adjust the arrangement sequence of subunits in the supercell along y axis shown in Table 1, which is the 2 nd case marked as "C2". The corresponding model is shown in Fig. 9(b).

The phase shift of each adjusted subunit along y axis is shown in Fig. 9(c). It can be found that the phase shift of the subunits for C2 still keeps linear changing due to the wave periodicity, but the phase shift of the supercell is doubled. This results in that the for the supercells in which the number of subunits is even, the phase distribution will have some small discontinuities after changing the arrangement sequence of the subunits in the above way, but the effect can be ignored when the phase resolution is much greater than .

By the above three cases, we show that adjusting the arrangement sequence of subunits in the supercell can suppress the 1 st order propagation diffraction. Based on this, we further show how to improve the supercell of the quasi-omnidirectional LGEM with a small thickness while maintaining its high-efficiency absorption. The thickness of LGEM with J subunits (J=13 in Fig. 9(a)), which has phase gradient of , is determined by the maximum subunit length. The subunit length consists of two parts according to Eq. ( 5). The first part is the length of , which produces linear gradient distribution of phase. The second part is the additional fixed length of , which is mainly to enhance the absorption of multiple reflections. Here, we aim to reduce the second part of the length. For reducing the first part, the corresponding detail can be found in Appendix E.

For the improved supercell, the subunit number is increased by times, while keeping subunit width, slit width, and additional fixed length unchanged. The subunit lengths are redesigned to keep a linear gradient phase distribution according to Eq. (5).

For example, the number of subunits in Fig. 9(a) is increased to , which is the 4 th case marked C4, as shown in Fig. 9(e). The subunit number is increased by times, but the phase gradient is reduced as . According to Eq. ( 17),

2 p 0 2k g > 0 2 j J l 0 h â =13 J =25 J ¢ ˆ2 a » 4 0
=1.05k g the 1 st order diffraction will dominate for the negative incident angle in the whole range from 343 Hz to 1000 Hz. This results in the low-efficiency absorption in this frequency range, the absorption coefficient of C4 is shown in Fig. 9(f).

Further, by adjusting the arrangement sequence of the above subunits, the phase shift of the supercell can be increased by times. The phase gradient and phase resolution of the above modified supercell can be rewritten as and z * , respectively. First, the modified phase gradient can satisfy the inequality of to suppress low-efficiency absorption induced by the 1 st order diffraction.

Through the inequality, we can get the relation of . Second, the modified phase resolution of z * can be less than to ensure good diffraction performance. We can get the increased times as . Therefore, when fulfil the relation of

, ( 32
)
the times of multiple reflections for 0 th order diffraction mode is increased by times, without introducing other propagation diffraction modes. In this way, the absorption effect of the modified supercell will be increased by times according to Eq. ( 27).

On the other hand, the additional fixed length of is reduced by times, the absorption effect from the additional length will be reduced by the same times.

Therefore, when the additional fixed length of is reduced by times, at the same time, the multiple reflections is increased by times in the above way, the modified supercell can still keep the same high-efficiency absorption as the unmodified one in C1.

For example, the arrangement sequence of the subunits in C4 is redesigned, which is the 5 rd case marked as "C5". When the subunit number is increased by times compared with C1, the phase gradient can be increased by times according to Eq. ( 32). In this way, the additional fixed length of can be reduced by times while maintaining the same high-efficiency absorption. The modified supercell structure with the additional fixed length of is shown in Fig. 9(g). The corresponding absorption coefficient is shown in Fig. 9(h). It can be found that in the frequency range from 343 Hz to 1000 Hz, the modified supercell can obtain the same high-efficiency absorption. Importantly, compared with the 1 st case, the modified supercell thickness in propagation direction is less than the wavelength in the whole frequency range from 343 Hz to 1000 Hz and only 0.57 times the maximum wavelength.

This confirms validity of the general approach, which can also be applied to the absorption of acoustic and electromagnetic wave metasurfaces. With these methods, we can further compress the thickness of the LGEM by a larger multiple to get a smaller or even deeper sub-wavelength thickness.

Experiments of flexural waves absorption by the LGEM

In order to further validate the vibrations absorption of the LGEM, we adopt the

LGEM designed in the section 4.1 to perform a set of experimental investigations. The

LGEM with phase gradient is composed of twelve different subunits, as shown in Fig. 10

ˆ2 a » ˆ4 b = 0 h ˆ2 a » 0 0 = 2 h h ¢ 0 =2.1k g i).
Then, the damping layer (butyl rubber with a thickness of 3mm) is bonded on the back of the strip-like aluminum plate, as shown in . Finally, the constrained layer (aluminum foil with a thickness of 0.15 mm) is bonded to the upper surface of the damping layer as shown in Fig. 10(b-iii).

An array of eight piezoelectric patches (20 mm × 20 mm × 0.3 mm) with a distance of from the LGEM is bonded on the surface of the plate as actuators. The LGEM and the piezoelectric array are obliquely placed. In this way, the upward propagating flexural wave excited by the piezoelectric patches is reflected by the upper boundary into the left area attached with blue-tack (Liu et al., 2017b). The reflected flexural wave will be fully absorbed in a large area of blue-tack. The boundary of blue-tack in the middle of the plate can be equivalent to the non-reflection boundary, as shown in Fig. 10(a), which provides a reliable condition for the later steady-state experiments.

The built-in generator in the PSV-400 scanning laser Doppler vibrometer generates a sinusoidal signal and transmit the signal to a power amplifier (HVPA05). The power amplifier amplifies the signal and drives the eight piezoelectric patches. For the six piezoelectric patches in the middle, the amplified voltage is two times of that for the two outer piezoelectric patches. It should be pointed out that for different excitation frequencies, one need to adjust the distance between piezoelectric patches to generate a flexural wave Gaussian beam according to the phased array theory (Zhu et al., 2014b;Cao et al., 2018c). For example, in the central frequency of 600 Hz, the distance between piezoelectric plates is 2.5 times length of the piezoelectric patch width. Before the experiment, we also use COMSOL Multiphysics 5.4 software ("Piezoelectric Six measurement points, marked in Fig. 10(a), are selected in the measurement area. These measurement points are divided into three groups: points 1 and 2, points 3 and 4, points 5 and 6. Their locations are marked in Fig. 10(a), which meets the farfield assumption [START_REF] Denis | Measurement and modelling of the reflection coefficient of an Acoustic Black Hole termination[END_REF]. Therefore, the measurement wave field in these points can be well approximated as the sum of incident and reflection propagating waves. The out-of-plane complex displacements of these wave fields are measured by the laser scanning head of PSV-400 in the measurement mode of the frequency domain, which is perpendicularly aimed at the plate surface. An ensemble average with 50

samples are used at every measurement point to ensure the signal quality.

First, for the points 1 and 2, we test the complex displacements, which can be expressed as and ,

respectively. The transfer function H12 of the total wave field in the two points can be obtained from the measured complex displacements of w 1 and w 2 . It can be expressed as , where the reflection coefficient . Then, we can calculate the reflection coefficient of by the transfer function H12, the distance s1, s2, and the wavenumber k. Further, the experimental absorption coefficient can be obtained as .

For the other two groups of measurement points, the absorption coefficient can also be obtained using the same method. Then, the average value of the three groups of measured absorption coefficients is the reliable reflection coefficient. These experimental results are shown in Fig. 11, which are in good agreement with both the theoretical and numerical simulation results. It shows that there is a high-efficiency absorption in the frequencies range approximately from 340 Hz to 1000 Hz. Some small For comparison, the experiment of flexural wave absorption by the LUEM is also conducted. The test piece is illustrated in Fig. 11, and its back is bonded with constrained damping layer (not shown in the figure). It can be seen that the absorption coefficient is obviously smaller than that of the LGEM. The results are in good agreement with both theoretical and numerical simulation results. The comparison confirms that the gradient phase greatly enhances the absorption of flexural waves. In order to show the absorption effect intuitively, the dynamic full wave field for the central frequency of 600 Hz is obtained in measurement mode of the time domain, as shown in Supplementary Video. It can be seen that the incident wave propagating in the host plate is not superposed by the reflected wave from the boundary to form a standing wave. It intuitively shows that the incident wave propagating to the boundary has been efficiently absorbed by LGEM.

( ) ( ) 1 2 1 2 1 I R i s s k i s s k w W e W e + × -+ × = × + × 2 2 2 I R is k is k w We W e × -× = + ( ) ( ) ( ) 1
é ù = = + + ë û ! ! R I / r W W = ! ( ) ( ) 1

Conclusions

We have expanded the conventional mode-coupling method to analytically obtain the full field reflection of each flexural wave diffraction mode, and revealed that the flexural wave absorption behavior is due to multireflection-enhanced absorption of 0 th or -1 th order diffractions. Further, we have introduced the concept of sub-wavelength

LGEM based on the 0 th order diffraction to achieve a broadband and quasiomnidirectional absorption of flexural waves. We reveal that low-efficiency absorption performance of the LGEM is mainly from the 1 st order diffraction, the impedance mismatch at the interface, and the degenerated gradient in the low-frequency. In addition, the general approach of adjusting the arrangement sequence of subunits in the supercell, which can eliminate low-efficiency absorption induced by the 1 st order diffraction, is proposed to further compress the thickness of the LGEM while maintaining its high-efficiency absorption. Finally, the experiment confirms the highefficiency absorption of the LGEM in the frequencies range approximately from 340 Hz to 1000 Hz. This research could open up possibilities in the family of elastic metamaterials and metasurfaces to realize broadband perfect absorption in lowfrequency by sub-wavelength structures. Eqs. (C4)-(C9), in turn, can be rewritten as In other words, the phase shift of adjacent subunits still keeps a larger value in high frequency, although the linear phase gradient has degenerated. This indicates that multireflection induced by the 0 th order diffraction in high frequency is not very dependent on linear phase shift.

D D D D K d K d d K d d d d d d d K K d d g é ù = + - + - ë û é - ù - ae ö - + + - ç ÷ ê ú + è ø ë û ( ) ( ) 30 1 2 0 1 2 0 2 3 0 1 0 0 1 2 2 = 2 c c c c c c c d K d K d K d g d K K K K K g ae ö - + + ç ÷ è ø + + + + 0 1 20 2 d d d + = 0 1 2 30 2 + 2 d d d d + = 1 2 1 2 1 eff (1 ) c c E g K d k u = + eff 0 1 2 = + + m m m m 0 0 0 m d r = 1 1 1 m d r = 2 2 2 m d r = 1 r 2 r ( ) 1/4 2 eff eff eff / k m D w = eff k eff D 0 0 0 0 1 2 2 2 2 0 0 0 0 0 0 0 0 3 3 3 3 0 0 0 0 0 0 0 0 1,1,1,1 , , , = , , , , , ,
ik ik k k T D k D k D k D k iD k iD k D k D k é ù ê ú - - ê ú ê ú - - - ê ú - - ê ú ë û eff eff eff eff 2 2 2 2 2 eff eff eff eff eff eff eff eff 3 3 3 3 eff eff eff eff eff eff eff eff 1,1,1,1 , , , = , , , , , ,
ik ik k k T D k D k D k D k iD k iD k D k D k é ù ê ú - - ê ú ê ú - - - ê ú - - ê ú ë û 1 T 2 T eff eff
÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ + = ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ è ø è ø è ø H q q q Τ H Τ Τ q q q H q q q ! ! ! ! " ! # # # # # # # ! ! # ! " # # # # ! ! ! ! ( ) [ ]
ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ + = ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ ç ÷ è ø è ø è ø H q q q Τ H Τ Τ q q q H q q q ! ! ! ! " ! # # # # # # # ! ! # ! " # # # # ! ! ! ! ( ) [ ] 4 2 2 2 2 , = , , , j n j n i i V V V V × - -× q

Fig. 1

 1 Fig. 1(a) shows the schematic of LGEM with predesigned geometries in the host

  imaginary wavenumbers . The real and pure imaginary wavenumbers represent propagating and evanescent flexural waves, respectively.

  FIG. 2 (a) The test set-up of Dynamical Mechanical Analysis (DMA). (b) The fitting curves for the

  . 3(b). The slope , bending moment M and shear force V can be expressed in term of the displacement as , , and , respectively, where the superscript indicates space differentiation with respect to the corresponding coordinate. The relationships between the state and coefficient vectors in Regions (I) and (II) can expressed respectively as , (10) where and are the transformation matrix between the state and the coefficient vectors. They are given in Appendix B.The boundary conditions at the interface, as shown in the illustration of Fig.3

  FIG.3 (a)and varying with frequency when the thicknesses of the constrained and

  . 4(b-ii) that the total number of the propagation diffraction modes has changed to 2 dues to the absence of the -1 th order mode.For intuitively displaying the above analytic results, the incident angles of -45 o and 45 o for the GEM and LGEM are chosen to show the full reflection fields. Based on the solved complex amplitudes of three diffraction modes including propagating and evanescent ones, the total full reflected wave fields in Region (P) can be solved by .(25)The analytical reflected wave fields for the GEM and LGEM are shown in and 4(b-iii), respectively. It can be seen that the reflection amplitudes for the GEM are 1 with the incident angles of -45 o and 45 o , while the reflection amplitudes for the LGEM are close to 0.8 and 0 with the incident angles of -45 o and 45 o , respectively. Furthermore, the corresponding simulation of full wave field are performed. It can be seen from Figs.

  FIG.4LGEM and GEM with surface phase gradient . (a-i) and (b-i) The corresponding

  FIG. 5 (a) The supercell of LGEM, which is composed of J subunits. (b) The supercell of LUEM

  FIG.6LGEM and GEM with surface phase gradient . (a-i) and (b-i) The corresponding

FIG. 7

 7 FIG. 7 (a) The model of the supercell of GEM. (c) The model of the supercell of UEM. The signals

From

  Fig. 8(a), we can see that for the large incident angles (>75°), the absorption coefficient sharply drops. In order to reveal this phenomenon, we set up an analytical model, as shown in the illustration of Fig. 8(b). The host plate and the waveguide subunits, which are marked as Regions (1) and (2), respectively. The dotted line is the interface between the two regions. The left and right ends of the model are the nonreflection boundaries. The coefficient vectors in Regions (1) and (2) can be expressed as and . At the interface of , according to the boundary conditions, we can get . Then, by deriving and simplifying, we get

  . (30), the reflection coefficients induced by impedance mismatch at the interface are calculated and shown by the hollow circle line of Fig. 8(b). The corresponding simulated results are described by the solid circle line and in very good agreement with the analytical results. For comparison, the reflection coefficients of the LGEM in Fig. 6(b-ii) are added, as shown in the black solid line of Fig. 8(b). It can be seen that reflection coefficients for the LGEM are almost the same as that for the

  FIG. 8 Performance evaluation of the LGEM. (a) The absorption coefficient varying with frequency

  (d) the 0 th order diffraction dominates in the frequency range from 200 Hz to 650 Hz, but in Fig. 8(a) the absorption coefficient sharply drops in the range from 200 Hz to 343 Hz. In order to explain this phenomenon, the change rate of phase shift for the supercell varying with the frequency can be shift for the supercell in the central frequency of 600 Hz.

  FIG. 9 (a) Phase shift of each subunit for five cases along y axis. (b) The supercell of the 1 st case is

  (a). The fabrication process is shown in Fig.10(b). First, the Wire Electrical Discharge Machining with a manufacturing precision of 0.05 mm is adopted to accurately fabricate the gradient strip-like aluminum plate by cutting a thin aluminum plate (2000 mm × 1000 mm × 3 mm). The back of the test piece is shown in Fig.10(b-

  FIG. 10 (a) The fabricated structure of the LGEM. The simulated excitation wave field is shown in

  FIG. 11 The absorption coefficients of the LGEM and the LUEM by experimental, analytical and

  On the other hand, it can be seen from Fig.8(c) that the value of increases with the frequency in the frequency range from 800 Hz to 1000 Hz. It means the linear phase gradient has also degenerated in the high frequency. To clarify the effect of the degenerated phase gradient in high frequency on the absorption, we study the absorption coefficients in the frequency from 1000 Hz to 2500 Hz with considering non-degenerated loss. The storage modulus and loss factors for different frequencies are fixed to that in 600 Hz. The corresponding absorption coefficients is shown in Fig.12(b). It can be found that there is a high-efficiency absorption in the whole frequency range, except for the low-efficiency absorption in the lower left corner. The lowefficiency absorption can be explained simply. The reason is that, in the frequency from 1000 Hz to 2500 Hz, three propagation diffraction modes, i.e., the -1 th , 0 th , and 1 st modes, will exist. In the lower left corner of Fig.12(b), the 1 st order diffraction dominates to result in the low-efficiency absorption. The high-efficiency absorption in the remaining region confirms that the degenerate phase gradient in high-frequency does not affect the 0 th order diffraction. It should be pointed out that the degenerate phase gradient in high-frequency is different from that in low-frequency. In the illustration of Fig.8(c), one can see the phase shift of the supercell for the frequency of 1000 Hz becomes bigger.

  d

Fig. 12

 12 Fig.12The absorption coefficient with considering non-degenerated loss in the frequency range (a)

Fig. 13

 13 Fig.13The reflection coefficients of incident wave varying with the thickness ratio and incident

=1 mm p

Acknowledgement

This project is supported by the National Natural Science Foundation of China (Grant No.11972296), the 111 Project (No. BP0719007), the Institute CARNOT (ICEEL), La Région Grand Est, the Fundamental Research Funds for the Central Universities (No. 310201901A005). The first author is grateful to the support of Innovation Foundation for Doctor Dissertation of Northwestern Polytechnical University (Grant No. CX201936) and China Scholarship Council (CSC No. 201806290176).

4 th case with the additional length of . (g) The modified supercell of the 5 th case with a short additional length of . (d), (f), and (h) The absorption coefficient varying with the frequency and incident angle for the 2 nd , 4 th , and 5 th cases, respectively.

Table 1 The arrangement sequence of the subunits along y axis Based on the above analysis, by changing the arrangement sequence of the subunits, the phase gradient can be increased by any positive integral multiple while keeping linear due to the wave periodicity. For example, the arrangement sequence of subunits is adjusted as the 3 rd case, which is marked as "C3". The corresponding arrangement sequence and phase shift of subunits along y axis are shown in Table 1 and Fig. 9(c), respectively. It can be found that the phase gradient of the C3 is increased by 3 times compared with that of C1, but the phase resolution z, which is

is reduced by 3 times. From the corresponding calculating results, we get that absorption coefficient of the C3 will be weakened in the high frequency of more than 700 Hz. This can be understood because good manipulation effect based phase requires high phase resolution [START_REF] Li | Theoryofmetascreen-basedacousticpassivephasedarray[END_REF]. The low phase resolution weakens the multi- 1)

(1 ) 12( 1)

12( 1) .

By applying the continuous boundary conditions of x-components of bending moment at the interface of x=0, and using the orthogonal relationship of the waveshapes and integrating with respect to y direction at the region , we can get ( ) ( ) ( ) . Therefore, the first part of the subunit length can also be reduced by changing the thickness of subunits =ϧ in the direction of z-axis, where ϧ is the thickness ratio. The performance of the subunit depends on the matching degree of the impedance at the interface between the host plate and the waveguide subunit. The better the impedance matching is, the smaller the reflection coefficient at the interface is. It will result in a better performance. The reflection coefficients of incident power varying with the thickness ratio and incident angle are shown in Fig. 13 with the white contour of 0.2. It can be found that changing the thickness ratio in this range of 0.4<ϧ<1, the reflection coefficients are less than 0.2. This means it is feasible to reduce the subunits length in this way. The aspects with more detail will be the scope of future works.