Broadband and high-numerical-aperture sharp focusing for waterborne sound with metagrating-based lens

Metalens with broadband and high-efficiency focusing functionality is desired in various underwater acoustic applications such as sonar and oceanography. Here we design and demonstrate a metagrating-based lens consisting of spatially sparse and wavelength-scale meta-atoms with optimized structures. With the help of grating diffraction analysis and intelligent optimization algorithm, the reflective metalens enables broadband and high-numerical-aperture focusing for waterborne sound over a 40 kHz-bandwidth for working frequency at 200 kHz. Full-wave numerical simulations unambiguously verify a sharp and high-efficiency focusing of sound wave intensity, with the full width at half maximum at the focal spot being smaller than 0.5λ and thus beating the Rayleigh–Abbe diffraction limit. Our work not only provides an intelligent design paradigm of high-performance metalens, but also presents a potential solution for the development of planar acoustic devices for high-resolution applications.

In this article, we utilize the principle of metagrating to realize broadband and sharp focusing of water-borne sound. In fact, high-efficiency manipulation of underwater acoustic waves has significant advantages in many application scenarios such as human health monitoring and underwater navigation, where electromagnetic waves are not adequate. For underwater acoustic applications such as sonar and oceanography, broadband and sharp focusing of sound signals is highly desired. It is even more attractive that such focusing effect can be achieved by a planar acoustic device with a simple structure. In the literature, the focusing functionality of acoustic wave was achieved by utilizing the gradient phase change along the metasurface [20,22]. To improve the focusing efficiency, some researchers used a hybrid metalens design [40], where the central region of the metalens is a gradient metasurface, and the surrounding regions are metagratings suitable for large deflection angles. More recently, people began to utilize metagratings to realize wave beam focusing for air-borne sound waves [41][42][43]. For example, Chiang et al designed a metalens consisting of scalable and simply-structured meta-atoms to focus ultrasonic waves at a single frequency [41].
Inspired by above works, in this article we combine the grating diffraction analysis with the intelligent optimization algorithm to design and construct a simply structured reflective metalens, which consists of spatially sparse and wavelength-scale meta-atoms. Broadband and high-numerical-aperture (NA) focusing for waterborne sound is achieved over a 40 kHz-bandwidth for working frequency at 200 kHz. Through full-wave numerical simulations, we unambiguously demonstrate a sharp and high-efficiency focusing of beam intensity, with corresponding full width at half maximum (FWHM) at the focal spot staying below 0.5λ/NA and beating the Rayleigh-Abbe diffraction limit.

A metagrating-based metalens
Let us begin with the right half part of the metalens, as shown in figure 1(a), which is composed of a series of metagratings. The full metalens can be obtained by combining the right half with its mirror (i.e., the left half). In the metalens, each meta-atom (the unit-cell of metagrating) is delicately designed so that it reflects the normally incident wave only along the pre-specified direction (i.e., toward the focal spot), as marked by the diffraction angle θ L (L = 1, 2, . . .) between the surface normal (dashed lines) and the reflected beam (red rays) by the Lth meta-atom. A zoomed view of the metagrating is shown schematically in the inset of figure 1(a), where each meta-atom is composed of an elliptical iron cylinder sitting on a rectangular base, backed by a sound hard plane. To achieve beam focusing, different meta-atoms must reflect the incident wave along different directions. Therefore, not only the period d of each local metagrating, but also the geometrical parameters of each meta-atom (including the major semi-axis a, minor semi-axis b, and rotation angle ϕ of the elliptical cylinder, as well as the width w and height h of the rectangular base) should vary from place to place across the metalens. In our study, the background material is water, with its mass density and acoustic velocity being ρ = 1000 kg m −3 and c 0 = 1490 m s −1 , respectively. The material of the elliptical cylinder and its base is iron, and the mass density, longitudinal and transverse velocities are ρ 1 = 7670 kg m −3 , c p = 6010 m s −1 , and c s = 3230 m s −1 , respectively.
According to the theory of grating diffraction, a plane wave normally incident on the grating will be reflected into several discrete diffraction orders. The direction of reflected wave is determined by the Bragg's condition [41], i.e., d sin α n = nλ 0 , where α n represents the diffraction angle for the nth-order reflected wave, and λ 0 = c 0 /f 0 is the wavelength in water at frequency f 0 . We demand that the −1st-order (n = −1) reflected wave (with α −1 = arcsin −λ 0 /d < 0) by all meta-atoms converge to the focal point (−Δx, F y ). To be more specific, the direction of −1st-order reflected wave by the Lth (L = 1, 2, 3, . . . ) meta-atom should satisfy the following conditions, where θ L = |α −1,L | = −α −1,L (so that θ L > 0) and d L represent the −1st-order reflection angle and grating period of the Lth metagrating, respectively. From equations (1) and (2), the reflection angle θ L can be expressed as a function of L, According to equation (3), after we specify the position of the focal point (−Δx, F y ) and the reflection angle θ 1 of the 1st meta-atom (L = 1), the reflection angles of other meta-atoms can be uniquely determined in a succeed way. To be more specific, when θ 1 = 16 • , the period of the 1st meta-atom is given by d 1 = λ 1 /sin θ 1 = 27.0 mm. Then the reflection angle θ L and grating period d L for other meta-atoms (L = 2, 3, . . .) can be subsequently determined, with corresponding curves as a function of L plotted in figure 1 Without loss of generality, we assume that the working frequency is f 0 = 200 kHz, so that the wavelength in water is λ 0 = 7.45 mm. We set F y = 15λ 0 at this frequency, then the horizontal position of the focal spot is given by Δx = 2.49λ 0 , as determined by equation (2). Thus, the absolute coordinates of the   focal point at f 0 are (−Δx, F y ) = (−18.5 mm, 111.8 mm). To obtain a high NA, the right half part of the metalens contains 27 meta-atoms, covering a wide reflection angle range from θ 1 = 16 • to θ 27 = 69.0 • . It is noted that in the scattering process of waterborne sound by meta-atoms made of iron cylinders, both longitudinal and transverse elastic wave modes in solid iron are involved, which is much more complicated than the corresponding scalar wave scattering problem (e.g., electromagnetic or acoustic wave scattering), and thus presents much richer physics. For example, multiple-scattering effects of sound wave between any two meta-atoms need to be rigorously considered. To resolve this challenge, we utilize the genetic algorithm (GA, one kind of intelligent optimization algorithm) [34] to smartly design each metagrating so that all meta-atoms work cooperatively to achieve a broadband and high-efficiency focusing.
In this section, the sound-hard plane is used as an exemplary configuration to demonstrate the focusing functionality of the reflective metalens. But the 'sound-hard plane' here is basically a 'boundary condition' to make sure that all incident wave energy is totally reflected without any transmission. We note that a similar sharp focusing effect can be achieved with other boundary conditions such as a sound-soft plane, and an example of such metalens can be found in appendix A.

Genetical-algorithm-assisted intelligent design
In applying the GA, we set the geometric parameters (a, b, ϕ, w, and h) of each meta-atom as the optimization variables, and the summation of the −1st-order reflectivity f j R −1 by each meta-atom over different frequency f j (here f j covers a frequency range of 40 kHz around f 0 ) as the optimization targets. With the help of GA, we can quickly get the optimized geometric parameters of all meta-atoms, as shown explicitly in table 1. The details of the efficiency calculation can be found in appendix B.
In the first five rows of table 1, the reflection angles are relatively small (θ L < 41 • ), and each set of geometrical parameters is suitable for only one meta-atom and works for only one reflection angle. In contrast, in the last five rows of table 1, the reflection angles are relatively large (θ L > 41 • ), so that each set of geometrical parameters is applicable to more than one meta-atom and suitable for more than one reflection angle. For example, the 6th and 7th meta-atoms (i.e., L = 6, 7) share the same geometrical parameters (as specified in the 6th row of table 1), and they reflect the incident wave along the directions of  According to the grating diffraction theory, the nth-order diffracted wave is a propagating mode when its diffraction angle |α n | < 90 • . For the 1st meta-atom, its −1st-order reflection angle (θ 1 = 16 • ) is along the desired direction, but all other propagating orders are not. Since the diffraction angle of the nth order is given by sin α n = nλ 0 /d, it is easy to know that there are seven propagating diffraction orders (n = 0, ±1, ±2, ±3) if α −1 = −16 • . By applying GA to the design of this meta-atom, we can achieve a high-efficiency reflection (R −1 > 95%) along −1st-order channel over a broad range of frequency f ∈ [189.2, 228.4] kHz, as shown clearly by the reflection spectra in figure 2(a). Frankly speaking, this close-to-unitary efficiency of R −1 is not easy to obtain, because we need to strongly and simultaneously suppress the reflection efficiencies of other six propagating orders including R 0 , R +1 , R ±2 , and R ±3 . The scattered pressure field at f 0 = 200 kHz is plotted in the inset of figure 2(a), where a nearly perfect wave-front pattern along the R −1 channel can be observed. In fact, this anomalous reflection effect is induced by the interaction between the incident pressure wave in water and stress field in the solid material. In figure 2(a) we also plot the first Cauchy stress invariant in the iron cylinder, where complicated sub-wavelength features of the stress field are observed.
For the 2nd meta-atom, its −1st-order diffraction angle is θ 2 = 25.9 • , and there are five propagating diffraction orders (n = 0, ±1, ±2) need to be controlled simultaneously. By utilizing the GA for its design, a high-efficiency anomalous reflection (R −1 > 95%) along the n = −1 channel is achieved in a frequency range of f ∈ [193.5, 228.3] kHz, while other diffraction orders are strongly suppressed. Numerically simulated scattered pressure field at f 0 = 200 kHz is shown in the inset of figure 2(b), confirming the nearly perfect wave-front pattern along the desired reflection direction.
In contrast, the remaining meta-atoms (L = 3, 4, . . . , 27) induce only three propagating orders (n = 0, ±1) at f 0 = 200 kHz because sin θ L > 1/2. Therefore, we only need to suppress two unwanted diffraction orders (n = 0 and 1). This is also the reason why we can use the same set of geometrical parameters for several neighboring meta-atoms and reflection angles, as shown in the last five rows in table 1, which greatly simplify the design of metalens.
As shown in table 1, only ten set of geometrical parameters are enough for the design of all 27 meta-atoms. To show the generality and flexibility of the design of the last 22 meta-atoms (L = 6, 7, . . . , 27) with large deflection angles (θ L > 43 • ), in figure 3 we show the reflection efficiency along the R −1 channel for different operating frequencies and different reflection angles. Different colors represent different values of R −1 , and corresponding contours show the distribution of R −1 in the (f , θ L ) parametric plane. Obviously, a highly efficient anomalous reflection R −1 > 90% is achieved in a broad frequency range f ∈ [190, 230] kHz for those meta-atoms. We note that the contributions of wave components from large deflection angles are important for the realization of a high NA and sharp focusing.

Focusing performance of the metalens
Next, we combine the 27 meta-atoms to construct the right half part of metalens. Through numerical simulations, we show in figure 4(a) the intensity of scattered pressure field |p| 2 when a plane wave is normally incident on the metalens. As we can observe, even with the right half part, we have already achieved a broadband beam focusing effect, covering a frequency range from 190 kHz to 230 kHz. For the incident wave at f 0 = 200 kHz, for example, the focus position is designed to be F = (−Δx, F y ) = (−18.5 mm, 111.8 mm), which is marked by the crossing point of two white lines in figure 4(a) and correspond to F y = 15λ 0 and Δx = 2.49λ 0 at 200 kHz. In accordance, the focal spot given by full-wave numerical simulation is located at F = (−15.0 mm, 117.0 mm), which shows a slight shift from F predicted by the diffraction theory analysis. This minor shift is expected because equations (1)-(3) are based on the assumption that each meta-atom is placed in an infinitely grating with a constant period d, but in a real metalens each meta-atom lies in a finite grating with the period d L varying from place to place across the metalens. A similar sharp focusing effect is also observed for other frequencies, as shown clearly in figure 4(a). We see that the position of the focal point (−Δx, F y ) vary with the working frequency f (or wavelength λ), which is consistent with equations (1)-(3).
Finally, we combine the right half metalen with its mirror copy to construct a full reflective metalens with Δx = 16.4 mm to achieve an optimized performance over the whole frequency range [190 kHz, 230 kHz]. The numerical aperture NA = 0.93 of the full-lens is defined by the maximum reflection angle, i.e., NA = sin θ L=27 = sin(69.0 • ), which is considerably large as compared with previous works. The intensity distribution of the scattered pressure field is shown in figure 4(b), which clearly shows a sharp focusing effect in the desired frequency range. Figures 4(c) and (d) show the cross-sectional intensity profiles along the x-axis at the focal spot and corresponding FWHMs. The maximum intensity at the focus is 133.7 Pa 2 at 200 kHz, with the amplitude of incident wave being 1 Pa. Thus, a huge energy concentrating effect is obtained at the focus, which is beneficial to energy harvesting applications. More importantly, the FWHM at the focal plane is 0.382λ at this frequency, which has substantially broken the Rayleigh-Abbe diffraction limit. It is worth noting that the FWHM always stays below 0.5λ in the frequency range [190 kHz, 230 kHz], which is never achieved for the focusing of waterborne sound to the best of our knowledge. In this frequency range, the minimum value of FWHM is only 0.356λ, which is achieved at f = 190 kHz.

Sub-diffraction focusing
It is well known that due to the diffraction effect, the resolution of conventional optics is limited by the Rayleigh-Abbe diffraction limit of 0.5λ/NA, where λ and NA are the working wavelength and NA, respectively. Recently, there has been a growing interest in developing far-field super-resolution optical devices which can achieve sub-diffraction hotspot of size smaller than the Rayleigh-Abbe diffraction limit [44,45]. Such sub-diffraction focusing effect is induced by the interference of coherent propagating waves in the far-field regime, and is related to the fact that band-limited functions are able to oscillate arbitrarily faster than the highest Fourier components they contain, a phenomenon called super-oscillations. Accordingly, people proposed a criterion (0.38λ/NA) for the hotspot size of superoscillation focusing.
In our work, the sub-diffraction focusing of waterborne sound is achieved at a focus length of 15λ 0 , thus in the far-field regime. Furthermore, according to the grating diffraction equation, d L sin α n = nλ 0 , we always have d L > λ 0 for all n = ±1th order diffracted wave, which is clearly shown in figure 1(b). From our design principle, all −1th order diffracted wave by meta-atoms on the right-half lens and all +1th order diffracted wave by meta-atoms on the left-half lens converge to the same focal point, as shown schematically in figure 1(a). All these ±1th order diffracted waves are propagating waves. Higher-order diffracted waves are evanescent waves which decay toward zero rapidly for a distance 15λ 0 from the metalens. Therefore, the sub-diffraction focusing effect is induced by the coherent interference of the ±1th order (propagating) diffracted wave by all meta-atoms in the metalens, with the underlying physics being attributed to the superoscillation effect.
To this end, in figure 4(d) we have plotted two dashed lines marking the Rayleigh-Abbe limit (0.5λ/NA = 0.538λ) and the superoscillation criterion (0.38λ/NA = 0.409λ), respectively, so that the graph is divided into three regimes. The regime above the Rayleigh-Abbe diffraction limit is defined as the subresolved regime, the regime below the 0.38λ/NA criterion is defined as the superoscillation regime, and the regime in between is defined as the superresolution regime [44].
At this stage, we also note that mathematically speaking, there is no fundamental limit on the hotspot size of the superoscillation effect. In principle, waves can be focused into a sub-diffraction superoscillatory hotspot of any shape and size beyond the Rayleigh-Abbe limit. However, there is no free-lunch. In super-oscillation fields, high resolution comes as a price of low efficiency, narrow field of view, and big sidebands [44,45].

Conclusions
In summary, the concept of metagrating is utilized to construct a reflective metalens for the focusing of underwater acoustic waves. Each meta-atom consisting of elliptical and rectangular iron cylinders is smartly designed to realize individual and separate control of each diffraction order of reflected wave. With the help of grating diffraction theory and intelligent optimization algorithms, the −1st order reflected wave by each meta-atom converges to the same focal spot, covering a wide range of reflection angle (from 16 • to 69 • ), and thus leading to a broadband and high-efficiency focusing effect. Full-wave numerical simulations verify a sharp and high-NA focusing for waterborne sound over a 40 kHz-bandwidth for working frequency at 200 kHz. Interestingly, the corresponding FWHM at the focal spot is less than 0.5λ/NA and breaks the Rayleigh-Abbe diffraction limit, which makes the metalens a feasible solution for the development of planar acoustic devices for high-resolution applications.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Reflective metalens with a sound-soft plane
In this work, we design a reflective metalens that can realize broadband and high-NA focusing for waterborne sounds. In section 2 of the main text, the sound-hard plane (implying a vanishing normal velocity field everywhere on the plane) is used as an exemplary configuration to demonstrate the focusing functionality of the reflective metalens. In other words, the 'sound-hard plane' is basically a 'boundary condition' to make sure that all incident wave energy is totally reflected without any transmission. We note that a similar sharp focusing effect can be achieved with other boundary conditions such as a sound-soft plane (implying a vanishing pressure field everywhere on the plane).
An example of such metalens is shown in figure 5. For a normally incident plane wave at f 0 = 200 kHz, the focal spot produced by the right-half metalens and the full metalens are plotted in figure 5(a). The FWHM of wave intensity at the focal plane is 0.370λ, exhibiting the desired sharp focusing functionality, as shown in figure 5(b). Meanwhile, the maximum intensity at the focus is 152.4 Pa 2 , with the amplitude of incident wave being 1 Pa, implying a huge energy concentrating effect. Of course, in the design process we have to re-optimize the geometrical parameters of all meta-atoms to achieve such focusing performance, and the optimized parameters are listed in table 2. We observe that the geometrical parameters optimized for a sound-soft plane are generally different from those optimized for a sound-hard plane.

Appendix B. Calculation of diffraction efficiency
We calculate the reflection efficiency through numerical simulations by using the acoustic-structure interaction module in COMSOL Multiphysics. As shown in figure 6, the left and right boundaries of the unit-cell are periodic boundary conditions, and the bottom boundary is set as a sound-hard plane. In addition, a perfectly matched layer is specified outside the upper boundary to absorb outgoing waves. For a normally incident plane wave, we calculate the reflection efficiency of each meta-atom by computing the integral of the pressure wave field on the upper boundary (marked as green line in figure 6) of the simulation domain. The pressure field of the normally incident plane wave is p in = b 0 e −ik 0 y , where k 0 = ω/c 0 is the wave-number in water. Thus, the energy flux along the y-direction is given by P in = 1 2 p in v * y,in dx = 1 2 −k 0 ωρ 0 |p in | 2 dx, where a negative sign '−' appears because the plane wave is incident towards the −y direction.
For the 0th order reflection wave, its wave vector is k r,0 = (0, k 0 ), and the corresponding pressure field and energy flux along the y-direction are given by p r,0 = b r,0 e ik 0 y , and P r,0 = 1 2 k 0 ωρ 0 |p r,0 | 2 dx, respectively. The reflection coefficient b r,0 can be calculated by numerical integration on the domain boundary with b r,0 = p * sc e ik 0 y dx, where p sc denotes the scattered pressure field, and '' * '' means complex conjugate. After obtaining b r,0 , we can calculate the energy flux P r,0 , as well as the reflection efficiency R 0 for the 0th order reflection wave through R 0 = P r,0 /P in .
For the ±1th order reflection waves, their wave vectors are given by k r,±1 = ±g, k 2 0 − g 2 , with g = 2π d . The pressure field and energy flux are, respectively, p r,±1 = b r,±1 e i(±g)x e i √ k 2 0 −g 2 y and P r,±1 = 1 2 √ k 2 0 −g 2 ωρ 0 |p r,±1 | 2 dx. The reflection coefficient b r,±1 can be calculated by numerical integration on the domain boundary with b r,±1 = p * sc e i(±g)x e i √ k 2 0 −g 2 y dx. After knowing b r,±1 , we can calculate the energy flux P r,±1 , as well as the corresponding reflection efficiencies through R ±1 = P r,±1 /P in .