When an acoustic wave with an incident angle θ i impinges on an arbitrary planar interface with two different media owning refractive indexes ( n i and n t ), the reflected angle θ r will be identical to θ i , and the refracted angle θ t will be determined by θ i , n i , and n t , which is the classical Snell’s law. However, when there exists a phase gradient at the interface, one can obtain the generalized Snell’s law: sin ( θ r ) − sin ( θ i ) = λ i 2 π d φ ( y ) d y , (1) sin ( θ t ) λ t − sin ( θ i ) λ i = 1 2 π d φ ( y ) d y , (2) where λ i and λ t are the wavelengths of the incident and transmitted wave, respectively. d φ ( y ) / d y is the phase gradient along with the planar interface. θ r and θ t can be expressed as: θ r = a r c s i n [ s i n ( θ i ) + ( λ i / 2 π ) ( d φ ( y ) / d y ) ] , (3) θ t = a r c s i n [ ( λ t / λ i ) s i n ( θ i ) + ( λ t / 2 π ) ( d φ ( y ) / d y ) ] . (4) Equations 3 and 4 imply that θ r and θ t are determined by three factors: incident angle, phase gradient, and wavelength. Thus, it is possible to manipulate the reflected and transmitted wave freely by designing a suitable phase shift φ ( y ) .

We consider a perforated unit cell composed of an air hole (the diameter being D) inside a square polymer (the width being a), as shown in Figure 1A. Figures 1B and C present the functional unit with nine periodically arrayed unit cells and BAAM, respectively. The blue medium indicates the water, and the gray and white medium are polymer and water. By arranging the functional unit with different diameters D, the metasurface can be constructed to steer the transmitted wave field. The material properties and geometric parameters of the BAAM are listed in Tables 1 and 2, respectively. In this article, the full-wave simulations are performed using commercial finite element software COMSOL Multiphysics with a preset Pressure Acoustic and Solid Mechanics module. Figure 2A illustrates the strip model with a functional unit. Continuous periodic boundary conditions are applied on the two sides of the strip model to simulate a periodically arranged structure. Perfectly matched layers (PMLs) are adopted at the upper end to absorb the reflection from the boundaries. The normally incident plane wave is employed at the bottom end of the strip model. To evaluate the phase shift and the transmission ratio, we set a probe line in the near area of the upper end, as shown in Figure 2A. The detected transmission ratio is defined as: T = | P o u t | | P i n | , (5) where p i n and p o u t are the amplitudes of the incident wave and the transmitted wave, respectively. Figure 2B shows the relationship between the phase shift and the diameter D. As D varies from 0 to 0.268 mm, the phase shift spans over a full 2π range, as shown with the black, and the transmittance exhibits a relatively high level, as shown with the color map. The number of air bubbles is an important variable for the functional unit. If the number of air bubbles is less than nine, it is impossible to achieve the phase shift covering an entire 2π range, and conversely, if the number of air bubbles is more than nine, the transmission ratio is affected, and the processing will be much more difficult.

Assuming that the acoustic wave is normally incident on the designed structure, Eq. 3 can be rewritten as: θ t = a r c s i n [ ( λ t / 2 π ) ( d φ ( y ) / d y ) ] . (6)

From Eq. 6, we can conclude that the transmitted angle varies with the phase shift, which indicates it is likely to steer the transmitted wave propagating along the direction as we want. Once we know the required phase shift and obtain the corresponding diameter D, the abnormal refraction will be achievable. This section performs the anomalous refraction effects of the proposed BAAM. Without loss of generality, we select three phase gradients: d φ = π 6 , 2 π 15 , π 9 . According to Eq. 6, we obtain the corresponding transmitted angles θ t = 20.9 o , 16.6 o , 13.8 o . For these three kinds of BAAM, Figures 3A–C illustrate the relationships between the phase shift and the diameter D. The blue dotted lines present the phase shift to manipulate the transmitted wave and determine the diameter D of each functional unit. The red dotted curves illustrate the diameter D corresponding to the phase shift. A Gaussian beam is normally applied to the bottom end of the metasurfaces as an initial incident at f = 700 kHz. As shown in Figures 3D–F, when the acoustic wave propagates through the metasurface, an apparent abnormal deflection of the acoustic wave can be observed in the transmitted field. The simulated results agree well with the theoretical results, which indicates the proposed metasurface shows an exceptional ability to manipulate the acoustic wave underwater. To verify the broadband feature of the proposed BAAM, we perform the full-wave manipulation at the frequencies of 694, 697, 700, 703, 706, and 709 kHz, respectively. Owing to limited space, we only demonstrate the simulated results of d φ = 2 π 15 with the transmitted angle θ t of 16.6°. As shown in Figures 4A–F, the proposed BAAM has the advantage of broadband in the abnormal refraction propagation.

In order to further verify the feasibility of the proposed metasurface to manipulate the acoustic wave, two different cases of self-bending beam propagation trajectories will be presented: the Bezier beam and the bottle beam.

In this section, we use a third-order Bezier curve to designate the acoustic propagation path of the transmitted wavefield. The Bezier curve is a polynomial function concerning the parameter t, which is expressed as: P ( t ) = ∑ i = 0 n P i n ! i ! ( n − i ) ! t i ( 1 − t ) n − i ( 0 ≤ t ≤ 1 ) , (7) where P(t) is the specific expression of a Bezier curve and n = 3 is the order. It can be found that once a series of points from P 0 to P n is selected, the trajectory of the curve will be determined. For a third-order Bezier curve, namely, a cubic Bezier curve, the expression will be obtained by defining four different points P 0 , P 1 , P 2 , and P 3 on the XY plane. Thus, the cubic Bezier curve can be expressed as a cubic polynomial: P ( t ) = ( 1 − t ) 3 P 0 + 3 t ( 1 − t ) 2 P 1 + 3 ( 1 − t ) t 2 P 2 + t 3 P 3                                                                 = ( 1 − t ) 3 ( 0 − 0.2311 ) + 3 t ( 1 − t ) 2 ( 0.1 − 0.2311 ) + 3 t 2 ( 0.25 0.1689 ) + t 3 ( 0.98 − 0.3311 ) . (8)

Figure 5A illustrates the process of obtaining the phase profile in the y-direction. By substituting different parameter t into Eq. 8, we obtain the designed trajectory (red curve) owning the shape of the Bezier curve with y = f ( x ) . Since the acoustic wave propagates along the direction perpendicular to the wavefront, we construct an envelope of tangential rays (black lines) of the trajectory, which will be utilized to build the wavefront (black dash line).

We present the tangential ray (red dash line) as an example, it intersects with the trajectory, wavefront, and y-axis at point A ( x , y ) , B ( u , v ) , and C ( 0 , ξ ) , respectively, where the distance of points B and C is δ , and the relationship between the two points satisfies the equations as follows: δ = − u cos ( π − θ ) = u cos ( θ ) , (9) ξ = v − [ − u t a n ( π − θ ) ] = v − u t a n ( θ ) . (10)

Since the acoustic rays are perpendicular to the wavefront, the relationship between the slope of points A and B can be expressed as: d x d z ⋅ d v d u = − 1. (11)

By simplifying Eq. 11, we can obtain: d u d v = − t a n ( θ ) , (12) where θ is the angle between the ray and the transverse direction. Lastly, with existing equations, the phase shift (blue line) in the y-direction can be expressed as: φ ( ξ ) = k δ . (13)

Combining Eqs 8–12, the relationship between the phase shift profile and the transmitted angle can be written as: d φ ( ξ ) d ξ = − k s i n ( θ ) = − k f ′ ( x ) 1 + f ′ ( x ) 2 . (14)

It can be found from Eq. 14 that once the desired beam path is determined, the phase shift to design the BAAM will be readily constructed. We use 31 functional units to build the BAAM to realize the Bezier beam effect. The simulated result of the acoustic pressure is provided in Figure 5B. The white dashed line represents the theoretical trajectory of the Bezier beam effect. It is easy to find that the simulated result is in good agreement with the theoretical value. Similar to the abnormal refracted, Figure 5C illustrates the desired phase distribution of 31 units obtained by Eq. 14. According to the relationship between the phase shift and the diameter D, the diameter D can be selected to design the BAAM consisting of 31 functional units. The distribution of the diameter D along the y-direction is shown in Figure 5D.

When the trajectory of the self-accelerating beam is a bottle trajectory, the transmitted acoustic wave propagates along a circular path, which is defined by: ( x − x 0 ) 2 + ( y − y 0 ) 2 = r 2 , (15) where ( x 0 , y 0 ) and r are the central and the radius of the circular trajectory, respectively. Substituting Eq. 15 into Eq. 14, we can obtain the phase profile of the bottle beam φ ( y ) along the y-direction: φ ( y ) = k t [ y − 2 r a r c t a n ( y / r ) ] , (16) where k t denotes the wavenumber of the transmitted wave. For the predesigned acoustic bottle beam, we set the center of the circle at ( x 0 , y 0 ) = ( 0.005 , 0 ) with the radius r being 5 mm. From Eq. 16, the phase shift of the circular trajectory can be further rewritten as   φ ( y ) = k t [ y − 0.01 arctan ( y / r ) ] .

Figure 6 (A) illustrates the schematic diagram of the circular trajectory described by y = 0.005 2 − ( x − 0.005 ) 2 . Since the functional units are arranged in a way of mirror-symmetrically around the x axis, the phase shift φ n = φ ′ n . Next, we assemble the metasurface with 30 functional arrayed units and test its performance with normal incident plane waves. The acoustic intensity of the transmitted wavefield is shown in Figure 6B. As can be clearly observed that the acoustic bottle beam propagates along the designated trajectory (white dashed circle) with nearly no sound pressure inside the bottle. The simulation results coincide well with the theoretical results. The phase shift can be obtained with Eq. 16, and the phase profile is shown in Figure 6C. Similarly, on the basis of the calculated phase shift from φ 1 to φ 30 , the diameter D can be selected to design the BAAM consisting of 30 functional units. The distribution of the diameter D along the y-direction is shown in Figure 6D.