# Refraction of Flexural Waves by Ultra-Broadband Achromatic Meta-Slab With Wavelength-Dependent Phase Shifts

Feng Liu, Zhichun Yang, Pengtao Shi, Yizhou Shen, Liyun Cao, Yanlong Xu

► To cite this version:

Feng Liu, Zhichun Yang, Pengtao Shi, Yizhou Shen, Liyun Cao, et al.. Refraction of Flexural Waves by Ultra-Broadband Achromatic Meta-Slab With Wavelength-Dependent Phase Shifts. Journal of Applied Mechanics, 2022, 89 (4), pp.041003. 10.1115/1.4053201. hal-05063018

HAL Id: hal-05063018

https://hal.science/hal-05063018v1

Submitted on 19 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Refraction of Flexural Waves by Ultra-broadband Achromatic Meta-slab with Wavelength-dependent Phase Shifts

Feng Liu<sup>a</sup>, Zhichun Yang<sup>a</sup>, Pengtao Shi<sup>a</sup>, Yizhou Shen<sup>a</sup>, Liyun Cao<sup>b</sup>, Yanlong Xu<sup>a*</sup>

<sup>a</sup> School of Aeronautics, Northwestern Polytechnical University, Xian 710072, China

<sup>b</sup> Universite de Lorraine, CNRS, Institut Jean Lamour, 54000 Nancy, France

E-mail: xuyanlong@nwpu.edu.cn (Yanlong Xu)

## Abstract

Great progress has been made in modulating flexural waves by elastic metasurfaces. Most of the proposed elastic metasurfaces suffer from chromatic aberration, limited in a narrow bandwidth around the designed frequency. In this paper, overcoming the chromatic aberration, an ultra-broadband achromatic meta-slab (UAM) with subunits of gradient thickness is proposed to realize the refraction angle unchanged with the incident frequency. Based on the phase compensation principle, wavelength-dependent phase shifts for the UAM that realize achromaticity are obtained. In order to verify the effectiveness of the theoretical design, the transmitted wavefields are solved according to the phased array theory, and the results correspond with those obtained by the finite element (FE) simulations and experiments, which show that the refraction angle is unchanged for incident frequencies from 2 kHz to 8 kHz. Besides, the UAM is extended into a periodic meta-slab, and multifrequency achromaticity is realized. Our designed meta-slabs overcome the chromatic aberration by simple configurations, which have significance in the applications of vibration control, vibrational energy harvesting, and health monitoring.

**Key words:** flexural wave, achromatic meta-slab, ultra-broadband, gradient thickness

## 1 Introduction

Chromatic aberration is a phenomenon that light rays of different wavelengths deflect at different angles or focus at different locations when they pass through a lens [1]. It significantly affects the functions of optical devices, such as imaging, communication, etc. In order to eliminate the chromatic aberration, the achromatic doublets, apochromatic, and superachromatic lenses are built by combining double, triple, and quadruple lenses [2-4] or
using the components with opposing dispersion properties [5]. Although the multifrequency achromaticity is realized, the combined lenses are complex and bulky, which are difficult and costly to apply.

Metasurfaces have an ultra-thin thickness and can abnormally modulate wave propagation by adjusting the phase shifts of subunits according to the generalized Snell's law (GSL) [6]. They provide a new method for eliminating chromatic aberration. By introducing an abrupt phase shift, various optical metasurfaces that realize achromaticity at multiple frequencies or in broadband have been demonstrated [7-13]. Recently, achromatic metasurfaces have been extended to acoustics. Tian et al. [14] designed a broadband metasurface with a gradient velocity and a gradient density, which exhibits a frequency-independent performance. Zhu et al. [15] proposed a specific gradient metasurface, whose phase gradient proportional to the wavenumber, to realize the dispersionless manipulation of reflected acoustic waves in broadband. Based on this, they developed a multifrequency achromatic acoustic metasurface [16]. Weng et al. [17] proposed a metasurface employed Helmholtz resonators, demonstrating its performance of achromatic acoustic focusing at multiple frequencies.

For elastic waves, various control forms are obtained by metasurfaces based on the GSL, including abnormal refraction [18-24], focusing [18,19,24-27], asymmetric transmission [28,29], total reflection [30], splitting [31], and other forms [23,32]. Typical reviews on metasurfaces refer to [33,34]. The metasurfaces above are designed for a single frequency. Once the metasurface configuration is established, the operating frequencies are limited in the band near the design one. In order to widen the frequency range, tunable metasurfaces are proposed [35,36]. However, achromatic metasurfaces for elastic wave modulations have not been paid attention. Flexural waves, a type of elastic waves in thin plates, have always been a concern because of the significant implications in wide engineering backgrounds [37-39]. It is also of great significance to eliminate the chromatic aberration of flexural waves. For instance, refracting flexural waves at the same angle or focusing flexural waves at the same location in broadband can effectively realize vibration control or vibrational energy harvesting.

We previously proposed a meta-slab with subunits of gradient thickness to realize abnormal refraction of flexural waves [21]. We have demonstrated that the phase shifts for the
subunits can be simply predicted by a formula related to the path length and wavelength. Based on this, in this paper, we design an ultra-broadband achromatic meta-slab with wavelength-dependent phase shifts according to the phase compensation principle, which can refract flexural waves at the same angle for different incident frequencies from 2 kHz to 8 kHz. The refraction angles obtained from the phased array theory, finite element (FE) simulations, and experiments correspond with the preset one, demonstrating the ultra-broadband achromaticity. Moreover, a multifrequency achromatic meta-slab (MAM) is realized by extending the UAM into a periodic structure. The multifrequency achromaticity is verified by the FE simulations, experiments, and analytical solutions solved based on the mode-coupling theory [40]. The proposed achromatic meta-slabs may push forward to apply to vibration control, energy harvesting, and health monitoring.

## 2 Theoretical design based on phase compensation principle

Consider that a beam of wave, having an arbitrary wavelength $\lambda$, is incident at an angle $\theta_i$ from point A, then passes through the meta-slab, and finally refracts out at an angle $\theta_t$ to point B, as shown in Fig. 1.(a). The total phase shift $\phi_{total}$ consists of three parts:

$$ \phi_{total} = \phi_{in} + \phi_{t} + \phi_{m} \qquad (1) $$

where $\phi_{in}$ and $\phi_{t}$ are the phases accumulated in the incident and transmitted regions, respectively, and $\phi_{m}$ is the phase shift imparted by the meta-slab. According to the geometrical relationship of the two infinitely close wave paths in Fig. 1(a), the difference of $\phi_{in} + \phi_{t}$, denoted by $d\phi_{in+t}$, is obtained as

$$ d\phi_{in+t} = \frac{2\pi}{\lambda} dy \cdot (\sin \theta_i - \sin \theta_t) \qquad (2) $$

In order to realize achromaticity, the beams of arbitrary wavelength should have the same wavefront at point B, that is, the total phase shift $\phi_{total}$ should be a constant independent of the wavelength. When $d\phi_{m}$ is equal to $-d\phi_{in+t}$, there is always $d\phi_{total} = 0$, indicating that the total phase shift is a constant independent of the wavelength, i.e., the achromaticity is realized.
Therefore, the phase shift gradient for the achromatic meta-slab is

$$d\phi_m = -d\phi_{in+t} = \frac{2\pi}{\lambda} dy \cdot (\sin\theta_t - \sin\theta_i) \quad (3)$$

Integrating Eq. (3) over y, we obtain the phase shift $\phi_m$, which is inversely proportional with the wavelength as follows:

$$\phi_m = \frac{2\pi}{\lambda} y(\sin\theta_t - \sin\theta_i) \quad (4)$$

However, the phase shift gradient for the conventional metasurfaces [18-20] and meta-slabs [21,31] is set as $d\phi_m/dy = 2\pi/p$. Then according to the generalized Snell's law: $\sin\theta_t - \sin\theta_i = (\lambda/2\pi)(d\phi_m/dy)$, the refraction angle is dependent on the operating wavelength and will change with the incident frequencies.

<figure><img src="image_2.png" /><figcaption>Fig.1. (a) Sketch for deriving phase shift for UAM. (b) Sketch of UAM embedded in host plate. (c) Enlarged view of UAM. (d) Phase shifts of subunits versus subunit thickness at 8 kHz.</figcaption></figure>

Based on Eq. (4), we propose an UAM that can deflect flexural waves with different wavelengths at the same angle. Figure 1(b) shows the sketch of the proposed UAM embedded in a host plate with a thickness of $h_0$. The UAM has the same neutral plane with the host plate.

Figure 1(c) shows the enlarged view of the UAM, which contains $m$ subunits with gradient thickness [$m=18$ in Figs. 1(b) and (c)], and the thickness of the $i$th subunit is denoted by
$h_i$. The subunits have the same width $a$ and are separated by the same slit with a width $s$.

The sum of the slit and subunit widths is denoted by $b$. The UAM length in the direction of wave propagation is denoted by $d$. The width of the UAM is denoted by $w$.

We have verified that the phase shift for the $i$th subunit can be directly expressed as [21]

$$ \phi_i = \frac{2\pi d}{\lambda_i} \qquad (5) $$

where $\lambda_i$ is the wavelength for the $i$th subunit. The phase shifts for the subunits with gradient thickness in Eq. (5) have the same form as in Eq. (4), which are both inversely proportional with the wavelength. Therefore, it is realistic to design an UAM using the subunits with gradient thicknesses. In our proposed UAM, the phase shift for the $i$th subunit according to Eq. (4) is

$$ \phi_i = \frac{2\pi}{\lambda}(\sin\theta_t - \sin\theta_i)ib \qquad (6) $$

According to Eq. (5) and (6), we obtain the difference of the phase shifts between the 1st subunit and the $i$th one as follows:

$$ \frac{2\pi}{\lambda_1}d - \frac{2\pi}{\lambda_i}d = \frac{2\pi}{\lambda}(\sin\theta_t - \sin\theta_i)(i-1)b \qquad (7) $$

where $\lambda_1$ is the wavelength in the 1st subunit. In addition, for flexural waves, there is $k_i = [\rho h_i (2\pi f)^2 / D]^{1/4}$, where $k_i$ is the wave number in the $i$th subunit; $f$ is the incident frequency; $\rho$ is the material density; $D$ is the flexural rigidity of the $i$th subunit. Combining the equation $\lambda_i = 2\pi/k_i$ with Eq. (7), we obtain the thickness of the $i$th subunit:

$$ h_i = \frac{d^2 h_0 h_1}{[d\sqrt{h_0} - \sqrt{h_1}(\sin\theta_t - \sin\theta_i)(i-1)b]^2} \qquad (8) $$

where $h_1$ is the thickness of the 1st subunit. Equation (8) indicates that the subunit thickness is independent of the incident frequency, that is, the ultra-broadband achromaticity is realized by using the meta-slab with subunits of gradient thickness.

It is well known that the flexural wave is a low-frequency approximation of the A0 Lamb wave, and there is a deviation between them when the product $fh$ of frequency $f$ and plate
thickness *h* is large. To demonstrate the reasonability of our design based on the thin plate theory, especially for the subunits with large thickness, we calculate the phase shifts of subunits at 8 kHz (the highest frequency investigated for UAM) for the subunit thickness ranging from 2 mm to 9 mm by using the thin plate theory and the Lamb wave theory, as shown in Fig. 1(d). It can be seen that the error between the phase shifts calculated by using the two theories increases gradually with the increasing of the subunit thickness. When the subunit thickness is less than 9 mm (*fh* < 72 kHz·mm), we notice that the error is less than 5%. This error does not influence wave modulations, which will be demonstrated in the following simulations and experiments.

Expect for the ultra-broadband achromatic deflection of flexural waves, the proposed UAM can also realize achromatic focusing. The broadband achromatic focusing performance has been demonstrated in our previous work [41] by using a shuttle-like meta-slab, whose achromatic principle is similar to the UAM.

## 3 Physical model based on phased array theory

To solve the full transmitted wavefield for the proposed UAM, a physical model for wave propagation is established according to the phased array theory. Since the subunit width is much less than the operating wavelength, the subunits can be regarded as a series of secondary point sources when the waves propagate out from them. Then the transmitted wavefield is a superposition of the radiation fields of these secondary point sources, as shown in Fig. 2.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Sketch of wave propagation based on phased array theory.</figcaption></figure>
The ith secondary point source is located at the interface between the ith subunit and the transmitted region, whose coordinate is $(x_i, y_i)$, and $y_i = (i-1)b + a/2$. The distance of an arbitrary point $(x, y)$ in the transmitted region from the ith secondary point source is denoted by $r_i : r_i = \sqrt{(x-x_i)^2+(y-y_i)^2}$. If there are m subunits, the full transmitted wavefield for an obliquely incident wave can be expressed as

$$w(x, y) = \sum_{i=1}^{m} e^{j\phi_i^0} (E_i H_0^{(2)}(k_i r_i) + F_i K_0(k_i r_i)) \quad (9)$$

where $\phi_i^0$ is the accumulated phase difference for the wavefront reaching the ith subunit relative to the reference point P; $E_i$ and $F_i$ are the transmitted ratios of the propagating and evanescent outgoing waves to the incident waves for the ith subunit, respectively; $H_0^{(2)}(k_i r_i)$ and $K_0(k_i r_i)$ are the second kind of Hankel and modified Bessel functions, respectively; $k_i$ is the wave vector of the transmitted waves.

For a beam impinging the meta-slab with an incident angle $\theta_i$, the wave path difference for the wavefront reaching the ith subunit relative to the reference point P is $l_i = y_i \sin\theta_i$, and thus $\phi_i^0$ is obtained as $\phi_i^0 = 2\pi l_i / \lambda$. The values of $E_i$ and $F_i$ can be obtained using the transfer matrix method, referred to our previous paper [21]. Note that $E_i$ and $F_i$ are complex coefficients and contain both the information of corresponding amplitudes and phase shifts for the ith subunit. So far, any parameter in Eq. (9) is known, and then the full transmitted wavefield can be obtained.

## 4 Analytical, numerical, and experimental demonstrations of ultra-broadband achromaticity

To demonstrate the ultra-broadband achromaticity, the commonly used 3D printing material polylactic acid (PLA), whose Young's modulus, density, and Poisson's ratio are 1.6916 Gpa, 1650 kg/m³, and 0.36, respectively, is adopted to print the UAM. The meta-
slab length *d*, subunit width *a*, slit width *s*, host plate thickness *h*₀, subunit number *m*
of the printed UAM, and the incident and refraction angles are chosen as 80 mm, 6 mm, 1 mm,
2 mm, 16, -5°, and 12°, respectively. The thicknesses for all subunits are obtained according
to Eq. (8) (see Table A1 in Appendix A for details). The achromaticity will be demonstrated by
full transmitted wavefields solved from the phased array theory introduced in Section 3, FE
simulations conducted in the Shell module of COMSOL Multiphysics v5.5, and experiments.

In the FE model, equivalent perfect matching layers are realized by setting their Young's modulus as 1.6716×(1+0.95i) Gpa, and their density and Poisson's ratio the same as those of PLA. A z-direction edge load with a magnitude of 1 is applied to generate plane flexural waves. The maximum mesh size is less than 1/15 of the minimum wavelength in the incident frequency range to ensure the accuracy of the FE simulation. The out-of-plane displacements are extracted to describe the full wavefields.

For experiments, the test piece is 3D-printed using PLA by the 3D printing technology.
The geometric parameters of the printed UAM are the same as those in the theoretical design.
The printed test piece is shown in Fig. 3(a), on which the sizes of the test piece are marked. An
array of six piezoelectric patches (10 mm × 8 mm × 0.3 mm) are bonded on the surface of the
test piece as excitations. The gap between adjacent piezoelectric patches is half of the patch
width. The voltage on each piezoelectric patch is the same to generate plane flexural waves.
Reflective membranes are pasted on the scanning areas of the piece to measure the signals
accurately. The incident angle of -5° is obtained through rotating the direction of the patch
array.
<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Test piece of UAM. (b) Experiment setup.</figcaption></figure>

Due to the limitation of the printing range of the 3D printer, the size of the test piece is small [see Fig. 3(a)], and the measurement aeras have to be set close to the boundaries of the test piece. Consequently, there are inevitable boundary reflections in the signals during the scanning in the time domain, which will introduce errors if the signals are transformed into the frequency domain by the fast Fourier transform. However, we can select the transient wavefields in the time domain before the boundary reflection happens to evade the refection impact. The 5-tune burst signal $V(t) = A_0[1 - \cos(2\pi f_c t/5)]\sin(2\pi f_c t)$ with the central frequency $f_c$ equal to the incident one is defined. It is generated by the function generator built in Polytec scanning laser Doppler vibrometer 400 (PSV-400) and amplified by a power amplifier (HVPA05). The vibrometer measures the out-of-plane displacements of the test piece.

The experiment setup is shown in Fig. 3(b).
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** (a) Polar directivity graphs of transmitted waves for frequencies with an interval 0.5 kHz from 2 kHz to 8 kHz. (b) Full wavefields, (c) transmitted wave fields, and (d) transient full wavefields obtained from FE simulations, phased array theory, and experiments, respectively, for frequencies 2 kHz, 4 kHz, 6 kHz and 8 kHz.</figcaption></figure>

The polar directivity of the transmitted waves in the frequency range from 2 kHz to 8 kHz are calculated by FE simulations, and the results for the frequencies with an interval 0.5 kHz
are shown in Fig. 4(a). They indicate that the refraction angle is unchanged with the incident frequency and identical to the preset one 12°. The full wavefields, transmitted wave fields, and transient full wavefields obtained from the FE simulations, phased array theory, and experiments, respectively, for incident frequencies 2 kHz, 4 kHz, 6 kHz and 8 kHz are shown in Figs. 4(b)-4(d). The refraction angles obtained from these three methods are in good agreement for each frequency and consistent with the preset one 12° signed by white arrows in all figures, which demonstrates the ultra-broadband achromaticity of the designed UAM.

In principle, the UAM can realize achromaticity regardless of the incident frequency. However, the achromatic performance, especially on the transmittance [15], is inevitably relevant to the frequency once the configuration is determined. We define the transmittance as $T = E_t/E_i$, where $E_i$ and $E_t$ are the incident and transmitted energies, respectively. The detailed calculation method of transmittance $T$ is given in Appendix B. We consider that $T > 0.6$ corresponds to the case that the UAM has a good achromatic performance. Since flexural waves are closely relevant to vibrations of plates, and their low-frequency propagation characteristics are more concerned in engineering, we only discuss the bandwidth when the incident frequency is less than 8 kHz.

The curves of transmittance versus the incident frequency for the UAMs with different number of subunits ($m=10, 14, 20$) but a fixed width $w$ [see Fig. 1(b)], i.e., 112 mm, are plotted in Fig. 5(a). The transmittances increase with the incident frequency and then tend to be flat, having consistent changing trends. Specially, the curves are almost overlapped, illustrating that the number of subunits has little influence on the transmittance when the width of the UAMs is fixed.

The curves of transmittance versus the incident frequency for the UAMs with different number of subunits ($m=10, 14, 20$) but a fixed subunit width $a$ (i.e., 6 mm) are shown in Fig. 5(b). Results show that the minimum effective frequency (corresponding to $T > 0.6$) increases with $m$ decreasing. In other words, when the subunit width is fixed, increasing the number of subunits is conductive to enhance the achromatic bandwidth. However, the number of subunits should not be too much, making the thicknesses of the subunits too thick and the interface impedances mismatch.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Transmittance versus incident frequency for UAMs with different number ($m = 10, 14, 20$) of subunits but (a) a fixed $w$ and (b) a fixed $a$.</figcaption></figure>

From Figs. 5(a) and 5(b), we can see that the transmittances are low for the frequencies less than 1 kHz. For these low frequencies, the wavelengths are relatively large, even larger than the width $w$ of the UAM and then the line source applied at the incident region of the FE model. Take 400 Hz for the red line in Fig. 5(b) as an example, the incident wavelength is 110.4 mm while the width $w$ of the UAM with ten subunits is 70 mm. In such a case the incident wave generated by the line source in the FE model behaves like a cylindrical wave, causing that the UAM cannot deflect the incident waves well. Therefore, the transmittance at this frequency is low. When the UAM width $w$ is increased, the transmittance at the same frequency will increase as well, as shown by the green and blue lines in Fig. 5(b).

# 5 Extension to multifrequency achromatic periodic meta-slab

The preceding designed UAM is an aperiodic structure. Extension of the UAM to a periodic structure can form a multifrequency achromatic meta-slab (MAM), which also breaks the single-frequency limitation of the existing metasurfaces or meta-slabs. In such a case, the UAM can be considered as a supercell of the MAM, and the supercell covers a full $2\pi$ phase shift range.

## 5.1 Multifrequency achromatic principle

The MAM, shown in Fig. 6(a), consists of periodical supercells (period length is $p$) arranged along the $y$-direction and has the same neutral plane with the host plate. The supercell contains $m$ [$m = 10$ in Fig. 6(a)] subunits with gradient thickness.
<figure><img src="image_7.png" /><figcaption>Fig. 6. (a) Sketch of MAM embedded in host plate. (b) Phase shift differences at frequencies $f_1$, $f_2$ and $f_3$.</figcaption></figure>

We consider flexural waves of a fundamental frequency $f_1$, and other frequencies, e.g. $f_q$ ($q=1, 2, 3, \dots$) incident to the MAM. Since the wavelength for a flexural wave is inversely proportional with the square root of the frequency, i.e., $\lambda \propto \sqrt{1/f}$, if

$$f_q = q^2 f_1 \qquad (10)$$

there is $\lambda^{(f_q)} = \lambda_1/q$, where $\lambda^{(f_q)}$ is the wavelength corresponding to $f_q$. Recalling Eq. (5), there is

$$\phi_i^{(f_q)} = q\phi_i^{(f_1)} \qquad (11)$$

Combining Eq. (5) with Eq. (11), we obtain

$$\frac{\lambda^{(f_q)} d\phi^{(f_q)}}{2\pi dy} = \frac{\lambda^{(f_1)} d\phi^{(f_1)}}{2\pi dy} \qquad (12)$$

According to the GSL, the refraction angle $\theta_i^{(f_q)}$ is

$$\begin{aligned} \theta_i^{(f_q)} &= \arcsin\left(\frac{\lambda^{(f_q)} d\phi^{(f_q)}}{2\pi dy} - \sin\theta_i\right) \\ &= \arcsin\left(\frac{\lambda^{(f_1)} d\phi^{(f_1)}}{2\pi dy} - \sin\theta_i\right) \\ &= \theta_i^{(f_1)} \end{aligned} \qquad (13)$$

We can see that the refraction angles $\theta_i^{(f_q)}$ for different frequencies $f_q$ are the same and equal to that for the frequency $f_1$.
Actually, the phase shift that the supercell spans for the frequency $f_q$ are $q$ times of $2\pi$.

According to Eq. (10), taking frequencies $f_1, f_2$, and $f_3$ as examples, the phase shift differences between adjacent subunits are shown in Fig. 6(b). The phases at the junctions of adjacent supercells are all matched (phases $2\pi$, $4\pi$ and $6\pi$ match the phase 0). Consequently, the designed meta-slab is achromatic at these frequencies. The MAM can be regarded as a periodic extension of the UAM. Specifically, the UAM has an increasing phase shift profile along the interface and do not constrained by the $2\pi$ phase range, while the MAM consists of supercells periodically along the interface, and the phase shift in the supercell is restricted by the full $2\pi$ phase range.

## 5.2 Analytical, numerical, and experimental demonstrations of multifrequency achromaticity

To verify the correctness of the above analyses, we derive the full transmitted wavefields by using the mode-coupling theory, from which, we obtain the refraction angle. The entire system shown in Fig. 6(a) can be divided into three regions: incident region (region I), meta-slab region (region M), and transmitted region (region T). The total displacement field in region I, formed by the superposition of the displacements of the incident wave and the reflected waves of all diffraction modes, can be expressed as

$$w^1(x, y) = A_i e^{-jk_{x0}y} e^{-jk_{x0}x} + \sum_n (A_n e^{jk_{xn}x} + B_n e^{\tilde{k}_{xn}x}) e^{-jk_{yn}y} \quad (14)$$

where $A_i$ is the amplitude of the incident wave, which is equal to 1 in this paper; $n$ represents the diffraction order; $A_n$ and $B_n$ are the reflection amplitudes corresponding to the propagating and evanescent flexural waves of the $n$th order diffraction mode, respectively; $j$ is the imaginary unit; $k_{x0} = \sqrt{k_0^2 - k_{y0}^2}$ and $k_{y0} = k_0 \sin \theta_i$ are the wave vectors of the incident wave in the $x$- and $y$- directions, respectively, with $k_0$ the wave vector of the incident wave and $\theta_i$ the incident angle; $k_{xn} = \sqrt{k_0^2 - k_{yn}^2}$ and $\tilde{k}_{xn} = \sqrt{k_0^2 + k_{yn}^2}$ are the wave vectors of the reflected propagating and evanescent flexural waves of the $n$th order diffraction mode in the $x$- direction, respectively, and $k_{yn}$ is the wave vector of the $n$th order diffraction waves in
the y- direction. According to the diffraction theory: $k_{yn} = k_0 \sin \theta_i + nG$, where $G=2\pi/p$ is the reciprocal lattice vector, with $p$ the periodicity of the supercell as shown in Fig. 6(a).

In region M, only the fundamental waveguide mode is considered due to the subunit width much smaller than the operating wavelength, and the displacement in the ith subunit is

$$w^M(i) = a_i e^{-jk_i x} + b_i e^{-k_i x} + c_i e^{jk_i x} + d_i e^{k_i x} \quad (15)$$

where $k_i$ is the wave number in the ith subunit; $a_i$, $b_i$ ($c_i$, $d_i$) are the coefficients of positive- and negative-going propagating (evanescent) flexural waves, respectively.

In region T, the displacement field, formed by the superposition of the displacements of the transmitted waves of all diffraction modes, can be expressed as

$$w^T(x, y) = \sum_n \left( C_n e^{-jk_{ny}x} + D_n e^{-\tilde{k}_{ny}x} \right) e^{-jk_{ny}y} \quad (16)$$

where $C_n$ and $D_n$ are the transmission amplitudes of the propagating and evanescent flexural waves of the nth order diffraction mode, respectively.

According to the interface continuity conditions at the interfaces of $x=0$ and $x=d$, we have

$$\left\{
\begin{array}{l}
w^I(0, y) = w^M(0, y) \\
\varphi^I(0, y) = \varphi^M(0, y) \\
\\
M^I(0, y) = \begin{cases} M^M(0, y), & |y - y_i| \le \frac{a}{2} \\ 0, & |y - y_i - b/2| \le \frac{b-a}{2} \end{cases} \\
\\
V^I(0, y) = \begin{cases} V^M(0, y), & |y - y_i| \le \frac{a}{2} \\ 0, & |y - y_i - b/2| \le \frac{b-a}{2} \end{cases}
\end{array}
\right.
\quad (17)$$

and
$$
\left\{
\begin{array}{l}
w^T(d, y) = w^M(d, y) \\
\varphi^T(d, y) = \varphi^M(d, y) \\
\\
M^T(d, y) = \begin{cases} M^M(d, y), & |y - y_i| \le \frac{a}{2} \\ 0, & |y - y_i - b/2| \le \frac{b-a}{2} \end{cases} \\
\\
V^T(d, y) = \begin{cases} V^M(d, y), & |y - y_i| \le \frac{a}{2} \\ 0, & |y - y_i - b/2| \le \frac{b-a}{2} \end{cases}
\end{array}
\right.
\qquad (18)
$$

where $\varphi^{\mathrm{I(T,M)}}$, $M^{\mathrm{I(T,M)}}$, and $Q^{\mathrm{I(T,M)}}$ are the slope, bending moment, and shear force of corresponding regions, respectively, which can be calculated by Eqs. (C1)-(C3) in reference [40].

Integrating Eqs. (17) and (18) with y and combining the orthogonal relation,

$$
\int_{0}^{p} e^{-jk_{1}y} e^{jk_{2}y} dy = 0 \quad (i_{1} \neq i_{2}) \qquad (19)
$$

we obtain the following matrix equation,

$$
\mathbf{X}[\mathbf{T}_r \quad \mathbf{T}_t \quad \mathbf{T}_m]^T = \mathbf{Y} \quad (20)
$$

where **X** is a matrix of size $[4m+4N] \times [4m+4N]$, and **Y** is a column vector with
$[4m+4N]$ rows, in which $m$ is the number of the subunits in a supercell and $N$ is the
number of total diffraction orders; $\mathbf{T}_r = [A_1 \ B_1 \ \cdots \ A_N \ B_N]$ is the vector of the reflection
coefficients of all diffraction modes in region I; $\mathbf{T}_t = [C_1 \ D_1 \ \cdots \ C_N \ D_N]$ is the vector of
the transmission coefficients of all diffraction modes in region T; $\mathbf{T}_M = [\mathbf{t}_1 \ \cdots \ \mathbf{t}_i \ \cdots \ \mathbf{t}_m]$
with $\mathbf{t}_i = [a_i \ b_i \ c_i \ d_i]$ is the vector of the coefficients of the fundamental waveguide mode
in Region M. The detailed derivation process is shown in Appendix C. All the coefficients of
$A_n, B_n, C_n$ and $D_n$ are solved from Eq. (20), and the displacement field in region T is
obtained according to Eq. (16). Then, we can obtain the analytical refraction angle from the
wavefield distribution.

To illustrate the multifrequency achromaticity, a MAM made of PLA is printed. The
fundamental incident frequency $f_1$ is randomly selected as 1 kHz, then $f_2$ and $f_3$ are 4 kHz and 9 kHz, respectively. The thickness of the host plate $h_0$ is 2 mm, and the corresponding operating wavelength at 1 kHz is 70 mm. The supercell length $d$ is 140 mm, which is twice the wavelength. The supercell contains ten subunits, and each subunit thickness is shown in Table A2 in Appendix A. The widths of the subunit and slit are randomly selected as 10 mm and 1 mm, respectively, and thus the theoretically predicted refraction angle is 39.5° according to Eq. (13).

Figures 7(a)-7(c) show the transmitted wavefields obtained from the FE simulations, mode-coupling theory, and experiments (the test piece of MAM is shown in Fig. D1 in Appendix D), respectively, for the perpendicular incidence with the incident frequencies of 1 kHz, 4 kHz, and 9 kHz, and the theoretically predicted refraction angle 39.5° is signed by white arrows in Figs. 7(a)-7(c) for comparison. The obtained wavefields are in good agreement, and the refraction angles read from them correspond with the theoretical one for each frequency. Besides, the polar directivity graphs for transmitted waves obtained by the experiments are given in Fig. 7(d), which can be seen that the refraction angles at these frequencies are unchanged, and the refraction angle 40.2° is almost equal to the theoretical one 39.5°. All results verify the multifrequency achromaticity of the designed MAM.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Transmitted wavefields at 1 kHz, 4 kHz and 9 kHz obtained from (a) FE simulations, (b) mode-coupling theory, and (c) experiments. (d) Polar directivity graphs for transmitted waves obtained</figcaption></figure>
by the experiments.

## 6 Conclusions

We analytically, numerically, and experimentally investigated the ultra-broadband achromaticity of the meta-slab with subunits of gradient thicknesses for refracting flexural waves. Based on the phase compensation principle, an UAM is designed, in which the phase shifts are wavelength-dependent, resulting in the subunit thickness wavelength-independent when the thickness of the 1st subunit is set to be the same as that of the host plate. The transmitted wavefields obtained from the phased array theory, FE simulations, and experiments correspond with each other, showing that the designed UAM deflects flexural waves at the same angle from 2 kHz to 8 kHz. When the subunit width is fixed, increasing the number of the subunits is conductive to reduce the minimum effective frequency and to enhance the achromatic bandwidth. Besides, extending the UAM into a periodic meta-slab can realize achromaticity at multiple frequencies, which is verified by the FE simulations, analytical solutions solved from the mode-coupling theory, and experiments. Our designs abstain from chromatic aberration by using simple configurations for the meta-slabs, potentially applying in vibration control, vibrational energy harvesting, and health monitoring.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant Nos. 12072276 and 11972296), and the 111 project of China (No. BP0719007).

## Appendix A Subunit thicknesses of meta-slabs

<table><caption>Table A1. Thicknesses of subunit in UAM (unit: mm)</caption><tr><th><i>h</i><sub>1</sub></th><th><i>h</i><sub>2</sub></th><th><i>h</i><sub>3</sub></th><th><i>h</i><sub>4</sub></th><th><i>h</i><sub>5</sub></th><th><i>h</i><sub>6</sub></th><th><i>h</i><sub>7</sub></th><th><i>h</i><sub>8</sub></th></tr><tr><td>2</td><td>2.11</td><td>2.22</td><td>2.35</td><td>2.49</td><td>2.64</td><td>2.80</td><td>2.98</td></tr><tr><th><i>h</i><sub>9</sub></th><th><i>h</i><sub>10</sub></th><th><i>h</i><sub>11</sub></th><th><i>h</i><sub>12</sub></th><th><i>h</i><sub>13</sub></th><th><i>h</i><sub>14</sub></th><th><i>h</i><sub>15</sub></th><th><i>h</i><sub>16</sub></th></tr><tr><td>3.18</td><td>3.39</td><td>3.63</td><td>3.90</td><td>4.20</td><td>4.53</td><td>4.91</td><td>5.33</td></tr></table>

<caption>Table A2. Thicknesses of subunit in MAM (unit: mm)</caption>
<table><thead><tr><th>h<sub>1</sub></th><th>h<sub>2</sub></th><th>h<sub>3</sub></th><th>h<sub>4</sub></th><th>h<sub>5</sub></th><th>h<sub>6</sub></th><th>h<sub>7</sub></th><th>h<sub>8</sub></th><th>h<sub>9</sub></th><th>h<sub>10</sub></th></tr></thead><tbody><tr><td>2</td><td>2.22</td><td>2.47</td><td>2.77</td><td>3.13</td><td>3.56</td><td>4.08</td><td>4.73</td><td>5.56</td><td>6.61</td></tr></tbody></table>

## Appendix B Detailed calculation method of transmittance

The transmittance is defined as the average energy ratio of the incident and transmitted waves. According to the fact that the energy is proportional to the square of the displacement amplitude, we can extract the average displacement amplitudes of the incident and transmitted waves (denoted by $w_i$ and $w_t$, respectively) to calculate the transmittance by using $T = (w_t/w_i)^2$. In detail, a strip with a span equal to one wavelength $\lambda$ and with an inclination angle equal to the refraction angle $\theta_i$ (the strip perpendicular to the wave vector) is defined in the transmitted region of the UAM. The average displacement amplitude calculated in this strip is $w_t$. Since the waves in the incident region of the UAM are the superposition of the incident and reflected waves, we instead define another strip in the incident region of the corresponding uniform plate, which has the span $\lambda$ and an inclination angle equal to the incident angle $\theta_i$, and calculate the average displacement amplitude in this trip as $w_t$.

## Appendix C Detailed derivation process of Eq. (20)

Integrating the first and second equations of Eq. (17) with $y$ within $|y-y_i| \le \frac{a}{2}$, we have

$$ 
\begin{aligned}
\mathbf{O}_1 + [\mathbf{P}_1^1 \quad \dots \quad \mathbf{P}_1^n \quad \dots \quad \mathbf{P}_1^N]^T \mathbf{T}_r &= \mathbf{Q}_1 \mathbf{T}_M \\
\mathbf{O}_2 + [\mathbf{P}_2^1 \quad \dots \quad \mathbf{P}_2^n \quad \dots \quad \mathbf{P}_2^N]^T \mathbf{T}_r &= \mathbf{Q}_2 \mathbf{T}_M
\end{aligned}
\quad (C1) $$

where

$$ \mathbf{O}_1 = \text{sinc}(k_{y0} \cdot a/2) \left[ e^{-jk_{y0}y_1} \quad \dots \quad e^{-jk_{y0}y_n} \quad \dots \quad e^{-jk_{y0}y_m} \right]^T, $$

$$ \mathbf{P}_1^n = \text{sinc}(k_{yn} \cdot a/2) \left[ \left( e^{-jk_{yn}y_1}, e^{-jk_{yn}y_2}, \dots, e^{-jk_{yn}y_n} \right)^T \right. \left. \dots \left( e^{-jk_{yn}y_1}, e^{-jk_{yn}y_2}, \dots, e^{-jk_{yn}y_n} \right)^T \right] $$

$$ \mathbf{Q}_1 = \text{diag} \left\{ \underbrace{[1 \quad 1 \quad 1 \quad 1]}_m, \dots, \underbrace{[1 \quad 1 \quad 1 \quad 1]}_m \right\} $$

$$ \mathbf{O}_2 = -jk_{x0} \text{sinc}(k_{y0} \cdot a/2) \left[ e^{-jk_{y0}y_1} \quad \dots \quad e^{-jk_{y0}y_n} \quad \dots \quad e^{-jk_{y0}y_m} \right]^T $$
$$
\begin{gather*}
\mathbf{P}_2^n = \text{sinc}(k_{yn} \cdot a/2) [\gamma_2^1 \ \dots \ \gamma_2^i \ \dots \ \gamma_2^m]^T \text{ with } \gamma_2^i = (jk_{xn} e^{-jk_{yn}y_i}, \hat{k}_{xn} e^{-jk_{yn}y_i})^T \\
\mathbf{Q}_2 = \text{diag} \left\{ \underbrace{[-jk_i - k_i \ jk_i \ k_i]}_{m}, \dots, \underbrace{[-jk_i - k_i \ jk_i \ k_i]}_{m} \right\}
\end{gather*}
$$

Integrating the third and fourth equations of Eq. (17) with $y$ within $0 \le y \le p$ and
combining Eq. (19), we have

$$
\begin{align}
\mathbf{O}_3 + \mathbf{P}_3 \mathbf{T}_r &= [\mathbf{Q}_3^1 \ \dots \ \mathbf{Q}_3^n \ \dots \ \mathbf{Q}_3^N]^T \mathbf{T}_M \tag{C2} \\
\mathbf{O}_4 + \mathbf{P}_4 \mathbf{T}_r &= [\mathbf{Q}_4^1 \ \dots \ \mathbf{Q}_4^n \ \dots \ \mathbf{Q}_4^N]^T \mathbf{T}_M
\end{align}
$$

where

$$
\mathbf{O}_3 = [D_0 p \cdot (-k_{x0}^2 - \nu k_{y0}^2) \quad \mathbf{0}]
$$

$$
\mathbf{P}_3 = D_0 p \cdot \operatorname{diag} \left\{ \left[ (-k_{x1}^2 - \nu k_{y1}^2), (\hat{k}_{x1}^2 - \nu k_{y1}^2) \right], \dots, \left[ (-k_{xN}^2 - \nu k_{yN}^2), (\hat{k}_{xN}^2 - \nu k_{yN}^2) \right] \right\}
$$

$$
\mathbf{Q}_3^n = p \sin c(k_{yn} \cdot a/2) [\mathbf{\zeta}_3^1 \ \dots \ \mathbf{\zeta}_3^i \ \dots \ \mathbf{\zeta}_3^m]^T \text{ with } {}_3\mathbf{\zeta}_i = D_i k_i^2 e^{jk_{yn} y_i} (-1, 1, -1, 1)^T
$$

$$
\mathbf{O}_4 = [D_0 p \cdot (jk_{x0}^3 + j(2-\nu)k_{x0}k_{y0}^2) \quad \mathbf{0}]
$$

$$
\mathbf{P}_4 = D_0 p \cdot \operatorname{diag} [\gamma_4^1 \ \dots \ \gamma_4^n \ \dots \ \gamma_4^N] \text{ with}
$$

$$
\gamma_4^n = \left[ (-jk_{xn}^3 - (2-\nu)jk_{xn}k_{yn}^2), (\hat{k}_{xn}^3 - (2-\nu)\hat{k}_{xn}k_{yn}^2) \right]
$$

$$
\mathbf{Q}_4^n = p \sin c(k_{yn} \cdot a/2) [\mathbf{\zeta}_4^1 \ \dots \ \mathbf{\zeta}_4^i \ \dots \ \mathbf{\zeta}_4^m]^T \text{ with } \mathbf{\zeta}_4^i = D_i k_i^3 e^{jk_{yn} y_i} (j, -1, -j, 1)^T
$$

Similarly, integrating the first and second equations of Eq. (18) with $y$ within
$|y-y_i| \le \frac{a}{2}$, we have

$$
\begin{align}
\begin{bmatrix} \mathbf{P}_5^1 & \cdots & \mathbf{P}_5^n & \cdots & \mathbf{P}_5^N \end{bmatrix} \mathbf{T}_r &= \mathbf{Q}_5 \mathbf{T}_M \tag{C3} \\
\begin{bmatrix} \mathbf{P}_6^1 & \cdots & \mathbf{P}_6^n & \cdots & \mathbf{P}_6^N \end{bmatrix} \mathbf{T}_r &= \mathbf{Q}_6 \mathbf{T}_M \nonumber
\end{align}
$$

where

$$
\mathbf{P}_5^n = \text{sinc}(k_{yn} \cdot a/2) [\gamma_5^1 \ \dots \ \gamma_5^i \ \dots \ \gamma_5^m]^T \text{ with } \gamma_5^i = (e^{-jk_{xn} \cdot d} e^{-jk_{yn} \cdot y_i}, e^{-\hat{k}_{xn} \cdot d} e^{-jk_{yn} \cdot y_i})^T
$$

$$
\mathbf{Q}_5 = \operatorname{diag} \left\{ \underbrace{\left[ e^{-jk_i \cdot d} \ e^{-k_i \cdot d} \ e^{jk_i \cdot d} \ e^{k_i \cdot d} \right]}_{m}, \dots, \underbrace{\left[ e^{-jk_i \cdot d} \ e^{-k_i \cdot d} \ e^{jk_i \cdot d} \ e^{k_i \cdot d} \right]}_{m} \right\}
$$

$$
\mathbf{P}_6^n = -\text{sinc}(k_{yn} \cdot a/2) [\gamma_6^1 \ \dots \ \gamma_6^i \ \dots \ \gamma_6^m]^T \text{ with } \gamma_6^i = (jk_{xn} e^{-jk_{xn} \cdot d} e^{-jk_{yn} \cdot y_i}, \hat{k}_{xn} e^{-\hat{k}_{xn} \cdot d} e^{-jk_{yn} \cdot y_i})^T
$$
$$
\mathbf{Q}_6 = k_i \cdot \operatorname{diag} \left\{ \underbrace{-j e^{-j k_i \cdot d} - e^{-k_i \cdot d} j e^{j k_i \cdot d} e^{k_i \cdot d}}_{m}, \dots, \underbrace{-j e^{-j k_i \cdot d} - e^{-k_i \cdot d} j e^{j k_i \cdot d} e^{k_i \cdot d}}_{m} \right\}
$$

Integrating the third and fourth equations of Eq. (18) with y within 0 ≤ y ≤ p and
combining Eq. (19), we have

$$
\begin{align}
\mathbf{P}_7 \mathbf{T}_r &= [\mathbf{Q}_7^1 \ \cdots \ \mathbf{Q}_7^n \ \cdots \ \mathbf{Q}_7^N]^T \mathbf{T}_M \nonumber \\
\mathbf{P}_8 \mathbf{T}_r &= [\mathbf{Q}_8^1 \ \cdots \ \mathbf{Q}_8^n \ \cdots \ \mathbf{Q}_8^N]^T \mathbf{T}_M \tag{C4}
\end{align}
$$

where

$$
\mathbf{P}_7 = D_0 p \cdot \operatorname{diag} \left[ \gamma_7^1 \ \cdots \ \gamma_7^n \ \cdots \ \gamma_7^N \right] \text{ with}
$$

$$
\gamma_7^n = \left[ e^{-jk_{xn} \cdot d} (-k_{xn}^2 - \nu k_{yn}^2), e^{-\hat{k}_{xn} \cdot d} (\hat{k}_{xn}^2 - \nu k_{yn}^2) \right]
$$

$$
\mathbf{Q}_7^n = p \sin c(k_{yn} \cdot a/2) [\boldsymbol{\zeta}_7^1 \ \cdots \ \boldsymbol{\zeta}_7^i \ \cdots \ \boldsymbol{\zeta}_7^m]^T \text{ with}
$$

$$
\boldsymbol{\zeta}_7^i = D_i k_i^2 e^{j k_{yn} y_i} (-e^{-j k_i \cdot d}, e^{-k_i \cdot d}, -e^{j k_i \cdot d}, e^{k_i \cdot d})^T
$$

$$
\mathbf{P}_8 = D_0 p \cdot \operatorname{diag} \left[ \gamma_8^1 \ \cdots \ \gamma_8^n \ \cdots \ \gamma_8^N \right] \text{ with}
$$

$$
\gamma_8^n = \left[ e^{-jk_{xn} \cdot d} (jk_{xn}^3 + (2-\nu)jk_{xn}k_{yn}^2), e^{-\hat{k}_{xn} \cdot d} (-\hat{k}_{xn}^3 + (2-\nu)\hat{k}_{xn}k_{yn}^2) \right]
$$

$$
\mathbf{Q}_8^n = p \sin c(k_{yn} \cdot a/2) [\boldsymbol{\zeta}_8^1 \ \cdots \ \boldsymbol{\zeta}_8^i \ \cdots \ \boldsymbol{\zeta}_8^m]^T \text{ with}
$$

$$
\boldsymbol{\zeta}_8^i = D_i k_i^3 e^{j k_{yn} y_i} (j e^{-j k_i \cdot d}, -e^{-k_i \cdot d}, -j e^{j k_i \cdot d}, e^{k_i \cdot d})^T
$$

When Eqs. (C1)-(C4) are simultaneous, we obtain

$$
X[T_r \quad T_t \quad T_m]^T = Y \quad (C5)
$$

where

$$
X = \begin{bmatrix} -\mathbf{P}_1 & -\mathbf{P}_2 & -\mathbf{P}_3 & -\mathbf{P}_4 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & -\mathbf{P}_5 & -\mathbf{P}_6 & -\mathbf{P}_7 & -\mathbf{P}_8 \\ \mathbf{Q}_1 & \mathbf{Q}_2 & \mathbf{Q}_3 & \mathbf{Q}_4 & \mathbf{Q}_5 & \mathbf{Q}_6 & \mathbf{Q}_7 & \mathbf{Q}_8 \end{bmatrix}^T, \quad Y = \begin{bmatrix} \mathbf{O}_1 & \mathbf{O}_2 & \mathbf{O}_3 & \mathbf{O}_4 & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix}^T
$$
## Appendix D Test piece of MAM

<figure><img src="image_9.png" /><figcaption>Fig. D1. Test piece of MAM</figcaption></figure>

## References

[1] E. Hecht, 1997, Optics, Addison Wesley 3.

[2] F.L. Pedrotti, L.S. Pedrotti, 1987, Introduction to Optics, Prentice-Hall.

[3] A. Miks, J. Novak, 2013, "Method for primary design of superachromats," Appl. Opt., **52**(28), pp. 6868-6876.

[4] G.I. Greisukh, E.G. Ezhov, S.A. Stepanov, 2006, "Diffractive–refractive hybrid corrector for achro- and apochromatic corrections of optical systems," Appl. Opt., **45**(24), pp. 6137-6141.

[5] C. Yuxin Wang, Ticknor, W. Yun, C. Jacobsen, 2003, "Achromatic Fresnel optics for wideband extreme-ultraviolet and X-ray imaging," Nature, **424**(6944), pp. 50-53.

[6] N. Yu, P. Genevet, M.A. Kats, F. Aieta, J.-P. Tetienne, F. Capasso, Z. Gaburro, 2011, "Light Propagation with Phase Discontinuities: Generalized Laws of Reflection and Refraction," Science, **334**(6054), pp. 333-337.

[7] G.H. Yuan, E.T. Rogers, N.I. Zheludev, 2017, "Achromatic super-oscillatory lenses with sub-wavelength focusing," Light Sci. Appl., **6**(9), p. e17036.

[8] S. Wang, P.C. Wu, V.C. Su, Y.C. Lai, C. Hung Chu, J.W. Chen, S.H. Lu, et al., 2017, "Broadband achromatic optical metasurface devices," Nat. Commun., **8**(1), p. 187.

[9] S. Wang, P.C. Wu, V.C. Su, Y.C. Lai, M.K. Chen, H.Y. Kuo, B.H. Chen, et al., 2018, "A broadband achromatic metalens in the visible," Nat. Nanotechnol., **13**(3), p. 227-232.
[10] M. Khorasaninejad, F. Aieta, P. Kanhaiya, M.A. Kats, P. Genevet, D. Rousso, F. Capasso, 2015, “Achromatic Metasurface Lens at Telecommunication Wavelengths,” Nano. Lett., **15**(8), pp. 5358-5362.

[11] O. Avayu, E. Almeida, Y. Prior, T. Ellenbogen, 2017, “Composite functional metasurfaces for multispectral achromatic optics,” Nat. Commun., **8**, p.14992.

[12] F. Aieta, M.A. Kats, P. Genevet, F. Capasso, 2015, “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” Science, **347**(6228), pp. 1342-1345.

[13] L. Huang, X. Chen, H. Muhlenbernd, G. Li, B. Bai, Q. Tan, G. Jin, T. Zentgraf, S. Zhang, 2012, “Dispersionless phase discontinuities for controlling light propagation,” Nano. Lett., **12**(11), pp. 5750-5755.

[14] Y. Tian, Q. Wei, Y. Cheng, Z. Xu, X. Liu, 2017, “Broadband manipulation of acoustic wavefronts by pentamode metasurface,” Appl. Phys. Lett., **107**(6), p., 221906.

[15] Y.F. Zhu, X.Y. Zou, R.Q. Li, X. Jiang, J. Tu, B. Liang, J.C. Cheng, 2015, “Dispersionless Manipulation of Reflected Acoustic Wavefront by Subwavelength Corrugated Surface,” Sci. Rep., **5**, p. 10966.

[16] Y.F. Zhu, X.D. Fan, B. Liang, J. Yang, J. Yang, L.L. Yin, J.C. Cheng, 2016, “Multi-frequency acoustic metasurface for extraordinary reflection and sound focusing,” AIP Adv., **6**(12), p. 121702.

[17] J.K. Weng, Y.F. Zhu, B. Liang, J. Yang, J.C. Cheng, 2020, “Wavelength-dependent multi-functional wavefront manipulation for reflected acoustic waves,” Appl. Phys. Express, **13**(9) p. 094003.

[18] J. Zhang, X.B. Zhang, F.L. Xu, X.Y. Ding, M.X. Deng, N. Hu, C.Z. Zhang, 2020, “Vibration control of flexural waves in thin plates by 3D-printed metasurfaces,” J. Sound Vib., **481**(19), p. 115460.

[19] J. Zhang, X. Su, Y. Liu, Y. Zhao, Y. Jing, N. Hu, 2019, “Metasurface constituted by thin composite beams to steer flexural waves in thin plates,” Int. J. Solids Struct., **162**, pp. 14-20.

[20] L. Cao, Z. Yang, Y. Xu, B. Assouar, 2018, “Deflecting flexural wave with high transmission by using pillared elastic metasurface,” Smart Mater. Struct., **27**(7), p. 075051.

[21] Y. Xu, L. Cao, Z. Yang, 2019, “Deflecting incident flexural waves by nonresonant single-
phase meta-slab with subunits of graded thicknesses,” J. Sound Vib., **454**, pp. 51-62.

[22] L. Cao, Z. Yang, Y. Xu, S.W. Fan, Y. Zhu, Z. Chen, B. Vincent, B. Assouar, 2020, “Disordered Elastic Metasurfaces,” Phys. Rev. Appl., **13**(1), p. 014054.

[23] H. Zhu, F. Semperlotti, 2016, “Anomalous Refraction of Acoustic Guided Waves in Solids with Geometrically Tapered Metasurfaces,” Phys. Rev. Lett., **117**(3), p. 034302.

[24] G. Su, Y. Zhang, Y. Liu, T. Wang, 2021, “Steering Flexural Waves by Amplitude-Shift Elastic Metasurfaces”, ASME J. Appl. Phys., **88**(5), p.051011.

[25] Z. Tian, L. Yu, 2019, “Elastic Phased Diffraction Gratings for Manipulation of Ultrasonic Guided Waves in Solids,” Phys. Rev. Appl., **11**(2), p. 024052.

[26] W. Wang, J. Iglesias, Y. Jin, B. Djafari-Rouhani, A. Khelif, 2021, “Experimental realization of a pillared metasurface for flexural wave focusing,” APL Mater., **9**(5), p. 051125.

[27] Y. Jin, W. Wang, A. Khelif, B. Djafari-Rouhani, 2021, “Elastic Metasurfaces for Deep and Robust Subwavelength Focusing and Imaging”, Phys. Rev. Appl. 15(2), p. 024005.

[28] B. Li, Y. Hu, J. Chen, G. Su, Y. Liu, M. Zhao, Z. Li, 2020, “Efficient asymmetric transmission of elastic waves in thin plates with lossless metasurfaces,” Phys. Rev. Appl., **14**(5), p. 054029.

[29] L. Cao, Y. Xu, B. Assouar, Z. Yang, 2018, “Asymmetric flexural wave transmission based on dual-layer elastic gradient metasurfaces,” Appl. Phys. Lett., **113**(18), p.183506.

[30] H. Zhu, T.F. Walsh, F. Semperlotti, 2018, “Total-internal-reflection elastic metasurfaces: Design and application to structural vibration isolation,” Appl. Phys. Lett., **113**(22), p. 221903.

[31] Y. Xu, L. Cao, P. Peng, X. Zhou, B. Assouar, Z. Yang, 2019, “Beam splitting of flexural waves with a coding meta-slab,” Appl. Phys. Express, **12**(9), p. 097002.

[32] Y. Liu, Z. Liang, F. Liu, O. Diba, A. Lamb, J. Li, 2017, “Source Illusion Devices for Flexural Lamb Waves Using Elastic Metasurfaces,” Phys. Rev. Lett. **119** (3), p. 034301.

[33] Y. Jin, Y. Pennec, B. Bonello, H. Honarvar, L. Dobrzynski, B. Djafari-Rouhani, M. I Hussein, 2021, “Physics of surface vibrational resonances: pillared phononic crystals, metamaterials, and metasurfaces”, Rep. Prog. Phys., **84**(8), p. 086502.

[34] Park, J., Lee, D., Rho, J., 2020, “Recent Advances in Non-Traditional Elastic Wave Manipulation by Macroscopic Artificial Structures,” Appl. Sci., 10(2), p. 547.
[35] S. Li, J. Xu, J. Tang, 2018, “Tunable modulation of refracted lamb wave front facilitated by adaptive elastic metasurfaces,” Appl. Phys. Lett., **112**(2), p. 021903.

[36] L. Ning, Y.-Z. Wang, Y.-S. Wang, 2020, “Active control of a black hole or concentrator for flexural waves in an elastic metamaterial plate,” Mech. Mater., **142**, p. 103300.

[37] X.-Q. Fang, J.-X. Liu, G.-Q. Nie, C. Hu, 2010, “Propagation of flexural waves and localized vibrations in the strip plate with a layer using Hamilton system,” Eur. J. Mech. A/Solids, **29**(2), pp. 152-157.

[38] H.F. Zhu, T.F. Walsh, F. Semperlotti, 2019, “Experimental study of vibration isolation in thin-walled structural assemblies with embedded total-internal-reflection metasurfaces,” J. Sound Vib., **456**, pp. 162-172.

[39] A. Zareei, A. Darabi, M.J. Leamy, M.-R. Alam, 2018, “Continuous profile flexural GRIN lens: Focusing and harvesting flexural waves,” Appl. Phys. Lett., **112**(2), p. 023901.

[40] L. Cao, Z. Yang, Y. Xu, S.-W. Fan, Y. Zhu, Z. Chen, Y. Li, B. Assouar, 2020, “Flexural wave absorption by lossy gradient elastic metasurface,” J. Mech. Phys. Solids, **143**, p. 104052.

[41] Y. Shen, Y. Xu, F. Liu, F. Wang, Z. Yang, 2021, “3D-printed meta-slab for focusing flexural waves in broadband,” Extreme Mech. Lett., **48**, p. 101410.