In the past decades, remarkable progress in metasurfaces has revolutionized wavefront shaping in a planar, miniaturized, and multifunctional way, which is enabled not by the phase accumulation during the propagation path but by the phase discontinuity at the interfacial boundary [1–3]. The abrupt phase discontinuity is encoded in elaborately designed subwavelength elements through various types of interaction with the incident light, such as electric and magnetic resonances [4,5], waveguide effect [6,7], and Pancharatnam–Berry phase [8]. The implementation of a continuous phase profile is through subwavelength discretization, and each element is only responsible for local phase modulation over its occupied subwavelength area. We name this type of wavefront shaping as local phase modulation (LPM), which is widely used in metasurface-based beam focusing [9–11], polarization control [12], holographic imaging [13,14], and so on.

Despite its remarkable success in wavefront control, LPM poses some limitations to metasurfaces. First, the phase response of each element is approximated by that of the element in locally periodic environment. This assumption fundamentally limits the application scenarios of LPM to slowly varying wavefronts. It further limits the material of metasurfaces to be high-index dielectrics or metals [15] as strong field localization is needed to minimize the coupling of adjacent elements in the aperiodic array. More importantly, LPM overlooks the impedance mismatch between the input and output beams, which inherently leads to increasingly strong parasitic scattering into undesired directions with the increase of the phase gradient [16]. Large-gradient phase profile is highly required in implementing tight focusing, beam multiplexing, or more complex wavefront shaping. In these cases, LPM-based metasurfaces become inefficient [6,17,18].

To circumvent these problems, metagratings [19] with superwavelength periods are proposed for near-perfect large-angle deflection by gathering energy into the first diffraction order (blazed diffraction) [20–25]. As constitutive elements of metasurfaces, metagratings with different periods and internal geometries are then combined for broad-angle beam scattering [26] and for metalenses with numerical aperture (NA) as high as 0.99 [27,28]. Each metagrating efficiently reroutes the beam following the ray path in geometric optics through efficient blazed diffraction.

However, the diffraction phase of metagratings is not studied yet. In fact, precise control of the diffraction phase is as important as maximization of the diffraction efficiency, as the incoming beam should be piecewise diffracted and constructively superimposed to build the desired wavefront. Furthermore, an arbitrary nonlinear-shaped large-gradient wavefront manipulation may not require piecewise blazed diffraction but other diffraction patterns. One should flexibly control the energy distribution among different orders and simultaneously control the diffraction phases so as to accurately reproduce a piece of arbitrary-shaped phase profile.

In this work, we propose a general mechanism for arbitrary large-gradient wavefront shaping by controlling a full set of diffraction parameters (direction, efficiency, and phase of each order) of superwavelength elements. The desired phase profile is cut into segments according to $2\pi$ zones, and the correlation between the phase segment and the desired diffraction parameters is revealed. The latter is then implemented by bilayer one-dimensional grating structures [29] with large design freedom. Different gratings are stitched into a metasurface to shape the wavefront in a desired way. As diffraction indicates the existence of lateral energy flow, we name this type of wavefront shaping mechanism as nonlocal diffraction modulation (nLDM).

As proof of the concept, we theoretically design and experimentally study high-NA cylindrical metalenses and 3D metalenses based on nLDM in the terahertz (THz) band. They show diffraction-limited focal spots, high efficiency, and subwavelength imaging resolution. As compared to LPM, nLDM serves as a better choice for large-gradient wavefront transformation. It brings additional merits of flexible choice of abundant dielectric materials and various processing methods, as strong field localization is not required. The method is general and can be extended for wavefront shaping at other frequencies.

Figure 1(a) shows an arbitrary wavefront shaping over a metasurface. The tangential field before and after the metasurface can be represented by $E_i(x,z=0_{-})$ and $E_o(x,z=0_{+})$. The ideal metasurface sitting at $z=0$ offers a transmission function $t(x)e^{j\phi (x)}$ to enable $E_o(x,z=0_{+})=t(x)e^{j\phi (x)}\times E_i(x,z=0_{-})$. $t(x)$ and $\phi (x)$ describe the amplitude and phase responses of the metasurface, whose thickness is overlooked here. In order to ensure high efficiency and to avoid the use of gain and lossy material, the metasurface is elaborately designed to implement the phase response $\phi (x)$, while keeping the amplitude as high as possible. Two design approaches of $\phi (x),$ the tradiational LPM and the nLDM, are schematically compared in Figs. 1(b) and 1(c). LPM discritizes the continuous phase with uniform subwavelength steps. Each step with a constant phase is realized by an element with proper topology. Different elements are stitched into the metasurface in the bottom of Fig. 1(b). At the positions where the phase variation is dramatic, such discritization is not accurate enough to reproduce the original curve. The phase delay of the subwavelength element is calculated in a periodic array, where the impedances of incident and transmitted beams are the same. But when the nonperiodic metasurface is constructed, the phase gradient will reshape the wavefront and change the transmission wave impedance. This impedance mismatch in the periodic array and nonperiodic array leads to parasitic diffraction of the energy into undesired directions, which become increasely severe in implementing large-gradient phase curves.

In contrast, nLDM divides the continuous phase curve into several segments, which can be described as $$\phi (x)=\{\begin{array}{cc}{\phi }_1(x)& x_1<x\le x_2\\ {\phi }_2(x)& x_2<x\le x_3\\ \mathrm{...}& \mathrm{...}\\ {\phi }_i(x)& x_i<x\le x_{i+1}\\ \mathrm{...}& \mathrm{...}\end{array}.$$(1)

There is no discritization error as compared to Fig. 1(b). The starting point and ending point ($x_i$ and $x_{i+1}$) of each segment can be chosen in different ways. As the phase profile is usually folded in $2\pi$, a segment can be simply chosen as a range from 0 to $2\pi$, as shown in Fig. 1(c). The length of the segments (${\mathrm{\Lambda }}_i$) varies in space and it will be slightly larger than the wavelengths in the large-gradient case. Please note that the choice of the segment can be flexible. If one segment is in the subwavelength scale, it can be merged into the latter one, as shown by the blue dashed line in Fig. 1(c). The starting and ending points can be arbitrary values not limited to 0 and $2\pi$. It is not straightforward to find a proper metagrating geometry to reproduce the exact phase segment. By considering the metagrating in an infinitely periodic configuration with ${\phi }_i(x)={\phi }_i(x+{\mathrm{\Lambda }}_i)$, the far field diffraction corresponds to the coefficients of the Fourier series expansion of the transmission function $e^{j{\phi }_i(x)}$ as [30] $$e^{j{\phi }_i(x)}=\sum _m{t}_m^d{e}^{j2m\pi x/{\mathrm{\Lambda }}_i},$$(2)where $t_m^d$ is the desired complex diffraction coefficient in the $m$-th order, which not only determines the diffraction efficiency as ${\eta }_m={|t_m^d|}^2$, but also defines the diffraction phase ${\psi }_m=\arg (t_m^d)$. Using Eq. (2), a piece of continuous phase profile in real space is mapped to $M$ diffraction efficiencies and $M$ diffraction phases, which can be used to find the target metagrating geometries. If the metagratings with desired diffraction behavior can be found, the phase profile can be implemented with high accuracy.

In order to meet so many diffraction requirements, we combine an optimization algorithm with electromagnetic simulation to find the metagrating quickly and automatically. We use the gradient descent optimization function fmincon in MATLAB and rigorously coupled wave analysis (RCWA) to approach the ideal diffraction by minimizing the following cost function: $$\mathrm{cost}=\sum _m|t_m^d-t_m^p(\mathit{g})|,$$(3)where $t_m^d$ and $t_m^p$ are the desired diffraction coefficients and the practical ones realized by the metagrating. $\mathit{g}=(g_1,g_2,g_3,\mathrm{...},g_N)$ represents the geometric parameters (widths of ridges and gaps) to be optimized. Since the optimization is not for the whole metasurface but its constitutive metagratings, it usually takes tens of seconds to find the optimized geometry. Equation (2) ensures that the ideal total diffraction efficiency is 100%, meaning that the beam should be totally transmitted and distributed among different orders. So high transmission is inherently considered during optimization, and the impedance mismatch issue is avoided. Coherent superposition of the segmented diffraction fields is realized by taking the diffraction phases into account.

Lastly, the metagrating should have enough degree of freedom so as to achieve different types of diffraction patterns. Here it is composed of bilayer one-dimensional ridges [Fig. 1(c)] made of polymer polylactic acid (PLA), whose refractive index is 1.57 at the target frequency 0.14 THz. Our previous study has proven the success of this type of metagrating for efficient large-angle diffraction [24,31], beam splitting [29], and asymmetric diffraction [32]. To ensure a planar profile of the metasurface, the thickness of each layer is not tuned. Throughout the manuscript, the ridges in the incident side are 2.0 mm thick, the slab is 2.1 mm thick, and the ridges in the transmission side are 3.1 mm thick if not mentioned.

As an example, we implement the red phase segment in Fig. 1(c) using LPM and nLDM at 0.14 THz. The length of this segment is $1.55\lambda$. The LPM design in Fig. 1(d) contains eight elements with fixed distance of 0.42 mm and spatially varied width. In the nLDM design, this phase profile is mapped to the diffraction coefficients of ($0.{88e}^{j0.5\pi }$, 0, $0.{35e}^{j0.8\pi }$) in the ($-1$, 0, 1) orders, which is achieved by the metagrating in Fig. 1(e). It contains only two ridges in each layer and has much larger feature size as compared to the LPM one. Simulations of the two designs are carried out in Lumerical FDTD Solutions. The field distributions in Figs. 1(d) and 1(e) show that the nLDM design can accurately reproduce the desired wavefront, while the LPM design fails.

As proofs of the concept, we design two high-NA cylindrical metalenses with one and four focal spots, respectively, using the nLDM method as shown in Figs. 2(a) and 2(b). The diameter of both metalenses is 82 mm ($37\lambda$). The focal length of metalens 1 is $f=15\mathrm{mm}$ ($7\lambda$), leading to a designed NA of 0.94. In metalens 2, four foci along the optical axis are located at $f_1=15\mathrm{mm}$, $f_2=25\mathrm{mm}$, $f_3=35\mathrm{mm}$, and $f_4=45\mathrm{mm}$. The ideal phase profile offered by the metalens can be calculated by backpropagation of point sources from the foci as $$\phi (x)=\arg \left(\sum _{n=1}^{n=N}{\xi }_n\times \frac{e^{j{k}_0{R}_n}}{4\pi R_n}\right),$$(4)where $k_0$ is the wave vector in vacuum, $R_n=\sqrt{f_n^2+x^2}$, and ${\xi }_n$ is the weight to balance the intensity ratio among different focal spots. $\xi =1$ for $N=1$ and ${\xi }_{1-4}=(0.3,0.6,0.8,1)$ for $N=4$ are used here. Figures 2(c) and 2(d) show half of the ideal phase profiles along the $x$ axis folded from 0 to $2\pi$. They are divided into segments whenever the phase reaches $2\pi$. Metalens 1 contains 14 segments in its half area and metalens 2 contains 11 segments. After Fourier transformation using Eq. (2), each phase segment is translated into a series of diffraction coefficients, whose amplitudes and phases are shown by the height and color of the histogram in Figs. 2(c) and 2(d). For metalens 1, all the segments require blazed diffraction with the same phase except for the first one. The largest blazed diffraction angle is 70°. The first phase segment is a quadratic shape, whose -1st order diffraction efficiency is 83% with the rest of the energy into the -2nd and 0th orders. The phase is also different from that of other metagratings. For metalens 2, the phase profile cannot be considered linear in several segments, leading to energy distribution in more than one diffraction order and variation of diffraction phases.

Next, using RCWA and the gradient descent optimization algorithm, we successfully find all the metagratings to reproduce the target diffraction patterns with $x$-polarized plane wave excitation (TM mode). Here we pick two segments in each metalens to show the optimization results, which are marked by the red phase segments in Figs. 2(c) and 2(d). The target and optimized diffraction efficiency (height) and phase (color) in each order is compared in Fig. 2(e), showing a good agreement. The phase profiles after the four optimized metagratings are plotted in Figs. 2(f)–2(i), which are consistent with the ideal phase segments quite well. The metagrating under each phase curve in Figs. 2(f)–2(i) shows that the relatively coarse geometry can reproduce the fast-varying phase profiles and even capture the deep subwavelength variation details as long as it exhibits the desired diffraction behavior. Details of the reproduction of phase profiles using the metagrating segments and the origin of the slight fitting error can be found in Appendix A.

The two metalenses are then 3D printed and experimentally tested. The metalenses are fabricated by a commercial 3D printer (Raise3D N2) using a fused-deposition modeling method. The diameter of the nozzle is 200 μm. The base plate and printing nozzle are set to 60°C and 215°C, respectively. The printing speed has a constant value of 60 mm/s. The metalens is piled up layer by layer along the y direction. Printing of each metalens takes 10 h approximately. The designed and 3D-printed metalenses are shown in Figs. 3(a) and 3(b), and the experiment setup in Fig. 3(c) is used to characterize the focusing performance. The source is a continuous-wave IMPATT diode at 0.14 THz. The beam is collimated by a Teflon convex lens with a focal length of 200 mm and a diameter of 4 inches. A half-wave plate is used to convert the polarization from the y to the x direction, so that the electric field is perpendicular to the ridges of cylindrical metalenses. To measure the field distribution after the metalens, a Schottky diode detector with an iris (diameter of 1 mm) is mounted on a three-dimensional translation stage and scanned with the step of 0.25 mm in the x, y, and z directions. For terahertz imaging, the imaging sample is mounted on the translation stage for scanning along the x and y directions in the focal plane. The detector is fixed at the optical axis after the sample for detection. A lock-in amplifier is connected to the source and the detector for modulation and lock-in detection.

For better comparison, simulations of the two metalenses are carried out in Lumerical FDTD Solutions with the Gaussian beam source polarized along the $x$ direction. The beam waist is 30 mm, which is consistent with that in the experiment. Perfectly matched layer boundary conditions are used in the x and z directions, and the periodic boundary condition is used in the y direction. The exit interface of the metalens is set as $z=0$. The simulated and measured focusing performances of the metalens 1 are shown in Figs. 4(a) and 4(b). The focal length is exactly 15.0 mm in simulation and 14.75 mm in measurement. The simulated FWHM of the $x$ component is 1.15 mm ($0.54\lambda$), equal to the diffraction-limited value. The FWHM becomes 1.55 mm ($0.72\lambda$) if the total intensity is considered. This is contributed by the strong longitudinal component away from the center due to the vector property of the field [33,34]. In experiment, the detector is responsive to the intensity of the total electric field, so the measured FWHM well agrees with the simulated one. Here we define the focusing efficiency as the ratio between the focused power and the input power, where the focused power is the transmitted power within the radius of 3 times of the FWHM in the focal plane [6]. Simulated focusing efficiency is 52% when $\mathrm{FWHM}=1.55\mathrm{mm}$ is considered. The measured focusing efficiency is 29.0%. Further study reveals that the discrepancy between the simulated and measured efficiency is likely caused by the fabrication error. The actual dimensions of the metalens show random deviation from the designed ones within $\pm 0.2\mathrm{mm}$. When the dimension error is considered for the metalens ridges and air gaps, the simulated focusing efficiency drops to 39% and the focal length remains fixed. So there is still room for improvement of the focusing efficiency in experiment by increasing the fabrication accuracy. Considering the simulated focusing efficiency above 40%, the focusing bandwidth covers 0.12–0.15 THz. The focal length is in linear proportion with the frequency, which is a signature of diffractive lenses.

In addition, the simulated focusing efficiency of 52% of the metalens 1 seems not a significant improvement if compared with some state-of-the-art metalenses in the optical region with similar NA designed using the LPM method [6,35,36]. However, they are all made of high-index dielectric on low-index substrate with very large aspect ratio, which is challenging at the terahertz band [10]. If we redesign the metalens 1 using the LPM method, the focusing efficiency is only 25.2%. So the 52% focusing efficiency does indicate the effectiveness of nLDM in boosting the performance of high-NA metalenses. In addition, the focusing efficiency can be improved to 59% and 71% in silicon-on-silicon metalenses and silicon-on-${\mathrm{SiO}}_2$ metalenses based on the nLDM method, as observed in Appendix B. Therefore, for any high-NA metalens with fixed material and fixed thickness, the nonlocal diffraction design in segments will always be more efficient than the local phase modulation of subwavelength elements. Another advantage is that the nLDM makes it possible to use very low-index dielectrics, such as plastic, glass, or photoresists, to build high-NA metalenses.

The simulated and measured field distributions through metalens 2 are shown in Figs. 4(d) and 4(e), which show excellent coincidence. Figure 4(f) plots the intensity along the optical axis, where four peaks can be identified. However, since the second and third peaks are very close, they cannot be distinguished through the measurement. Still the general trend and the relative intensity of the foci agree well with the simulation.

To further generalize the functionality from 2D focusing into full-space 3D focusing, the cylindrical metalens can simply be rotated around its center forming the concentric ring metalens. To ensure good quality of 3D tight focusing in the concentric ring metalens with linearly polarized excitation, it is necessary to have polarization-insensitive 2D metagratings, meaning the same diffraction patterns for both TE and TM polarization states. Since the metalens 1 is only optimized for TM polarization, we redesigned all the 14 segments in Fig. 2(c) by minimizing the cost function in both polarization states as $\text{cost}=\sum [|t_m^d-t_{m\_\mathrm{TM}}^p(\mathit{g})|+|t_m^d-t_{m\_\mathrm{TE}}^p(\mathit{g})|]$. As shown in Fig. 5(b), the optimized diffraction efficiencies in the dominate order ($-1$st order) for TE and TM polarizations are close to each other in each segment, and the diffraction phases of both polarizations are close to the target value $\pi$. The thicknesses of the bottom ridges, the slab, and the top ridges are 1.6 mm, 3.5 mm, and 1.6 mm, respectively. These polarization-insensitive metagratings are stitched and rotated into the metalens 3 in Fig. 5(a). This metalens is different from traditional Fresnel lenses realized by varying the thickness [37], as a full set of diffraction efficiencies and phases in each $2\pi$ zone are finely tuned. The focal spot is simulated and measured in Figs. 5(c) and 5(d). The polarization distribution of the focus is superimposed on the pseudo color of Fig. 5(c), showing linear polarization the same as the excited one.

To characterize the imaging performance of the metalens, we mount the imaging sample on a translation stage and scan it in the focal plane of the 3D metalens. A Schottky detector with a pin hole is fixed behind the sample to record the field distribution. The sample in Fig. 5(e) is a metallic sheet with etched “NK” letters whose linewidth is 2 mm. The detected intensity distribution can well resolve the letters. Next, we change the imaging sample into a resolution chart in Fig. 5(f). From left to right, the width of the bar is 2.0 mm, 1.8 mm, 1.4 mm, 0.9 mm, and 0.6 mm. The detected image shows a resolution of 1.4 mm, which is $0.65\lambda$. This agrees with the diffraction limit for $\mathrm{NA}=0.94$ considering the vectorial property.

We notice that each segment is designed and optimized in a periodic environment, while the whole device is aperiodic. This will inevitably bring some error. If we take the segmented metalens 1 as an initial solution and further numerically optimize the entire structure, the focusing efficiency is slightly boosted from 52% to 60%, indicating the rationality of the periodic assumption. This assumption now is in the superwavelength scale, therefore being better satisfied than the subwavelength periodic assumption in the LPM. For the case when the phase segment is drastically different from the neighboring ones, one can merge it to the neighboring one to form a larger segment. It will make the wavefront shaping more accurate. Nevertheless, the periodic assumption makes it possible to translate a continuous phase segment into a few discrete diffraction coefficients, which significantly reduce the information dimension, and simplify the optimization.

We know the LPM-based metasurface can be readily generated regardless of the device size once some basic subwavelength elements are found. In the nLDM, the optimization time will increase linearly with the scale of the metasurface, since each segment is optimized separately. However, the average optimization time for each segment is around 7 s. It is not a problem to deal with large-scale devices. In contrast, global optimization of the full metalens, though being a good solution to push the NA and to improve the efficiency [38–41], needs extensive computational resources and severely limits the scale of the devices. The proposed nLDM strategy is much more time-efficient in design and provides rich physical insights into the wavefront transformation.

Finally, the nonlocal diffraction engineering method is general enough to be applied at other frequency bands. When working at higher frequencies, it is possible to etch both sides of the glass plate or to build bi-layer photoresist structures similar to the one demonstrated here. Instead, one can simply choose high-index dielectric monolayer structures for nonlocal diffraction engineering. In Appendix B, we also give a monolayer metalens made of silicon with improved efficiency when designed using the nonlocal diffraction engineering, which can be scaled to optical frequencies and fabricated by conventional lithography techniques.

To summarize, we demonstrate that arbitrary large-gradient wavefront shaping can be accurately and effectively achieved from the nonlocal diffraction engineering of the superwavelength-scale metagratings. Different from the local phase modulation where each phase point corresponds to a subwavelength element, the proposed method here establishes correlation between the phase segments and the diffraction patterns, and further finds the proper metagrating geometry to fulfill the diffraction properties including diffraction angles, diffraction efficiencies, and diffraction phases. High-NA cylindrical metalenses composed of a series of metagratings are designed and experimentally tested, showing the effectiveness and flexibility of this method for arbitrary 2D wavefront shaping. The design is further extended to 3D metalenses by utilizing polarization-insensitive metagratings, whose application in high-resolution imaging is well characterized. These findings indicate that the nonlocal diffraction engineering offers a new perspective on beam transformation over metasurface devices, well suited for large-gradient cases.