By virtue of nonionization and good penetration for nonpolar materials, terahertz (THz) waves have been widely adopted in diverse applications, e.g., nondestructive testing^[1], biomedical imaging, and spectroscopy^[2,3]. High spatial resolution is pursued continuously in the prosperous THz imaging community. For continuous-wave full-field imaging, e.g., THz digital holography and THz ptychography, efforts have been made to mechanically or numerically expand the aperture of the array detectors in squeezed recording distance^[4–7]. For other imaging geometries using more sensitive cell detectors, their lateral resolution mainly depends on the diameter of the focusing spot of the illumination beam^[8–12], which is usually converged to $\sim 2\mathrm{mm}$ at minimum by off-axis parabolic mirrors or homogeneous transmissive convex lenses (TCLs)^[13,14]. The Abbe limit can be overcome by using composite lens made of highly resistive silicon and high-density polyethylene, so that the resolution is enhanced from $0.85\lambda$ to $0.35\lambda$ at 0.5 THz^[15]. The diameter of the incident beam is converged into $0.5\lambda$ at the focal plane by using a super-surface lens consisting of a stacked graphene strip within two layers of graphene^[16] or even $\lambda /150$ at 0.1 THz by utilizing THz surface plasmon polaritons^[17].

There is another category of components that not only shrinks the full width at half-maximum (FWHM) of the main lobe of the focal point but also extends its focal depth, which is critical when penetrating the thick barrier in many THz application scenarios. For instance, Bessel zone plate generates a THz beam with a focal depth of $50\lambda$ and an FWHM of $0.9\lambda$ at 0.6 THz^[18]. Various super-oscillation lenses (SOLs) have been fabricated at the THz band, among which phase-type and amplitude-type SOLs achieve a focal depth of $20\lambda$ and an FWHM of $1.2\lambda$ at 2.52 THz^[19], and a focal length of $25\lambda$ and an FWHM of $0.5\lambda$ at 0.1 THz^[20], respectively. However, the accompanying strong sidelobes affect resolution and fidelity. In addition, these components are usually fabricated using photolithography and laser ablation^[18–20]. With the rapid development of three-dimensional (3D) printing^[21], various diffractive elements, e.g., phase gratings^[22], spiral phase plates^[23], and vortex phase plates^[24], have been made using resin, acrylonitrile butadiene styrene plastic, and other polymers. Their transmittance, diffraction efficiency, and micrometer processing accuracy all satisfy the demand at the frequency band lower than 1 THz. Nevertheless, the research on using 3D printing to fabricate binary diffractive lenses (BDLs) with a subwavelength focal spot has not yet been reported. In this Letter, a BDL made up of concentric rings is proposed to realize a subwavelength focusing at 278.6 GHz. Its design, fabrication, characterization, and imaging applications are introduced in turn.

Assume that this diffractive element is illuminated by a monochromatic parallel wave at a wavelength of $\lambda$. The step height difference $\mathrm{\Delta }h$ between all concentric rings and their substrate satisfies a fixed phase delay of $\pi \text{:}\mathrm{\Delta }h=\frac{\lambda }{2(n-1)}$, where $n$ is the refractive index of the homogeneous material. Considering the circular symmetric structure of this BDL, the function of its complex transmittance can be set to $T(r)$, where $r$ is the distance to the center. Meanwhile, this structure is conducive to uniform beam shaping. According to the Rayleigh–Sommerfeld diffraction integral^[25], the intensity distribution of the optical field downstream from the BDL is also circularly symmetrical, $$\begin{gathered} I(r^{\prime })=|\frac{k}{iz}{\int }_0^{\rho }T(r)r\exp \left(ik\sqrt{{r^{\prime }}^2+r^2+z^2}\right) \\ {\times J_0\left(\frac{k{r}^{\prime }r}{\sqrt{{r^{\prime }}^2+r^2+z^2}}\right)\mathrm{d}r|}^2, \end{gathered}$$(1)where $r^{\prime }$ is the lateral distance to the optical axis, $k=2\pi /\lambda$ is the wavenumber, $z$ is the axial propagation distance, and $J_0$ is the zeroth order Bessel function.

According to Eq. (1), this diffractive element with proper structure is available to modulate the point spread function within the size of an Airy spot, thus providing a focusing capability of subwavelength magnitude^[26,27] and even manipulating the depth of focus^[26,28,29]. The radii of the rings are related with the optimization of the cost functions in simulated annealing algorithm, which is a heuristic random search process based on the solution strategy of Monte Carlo iteration^[28,30].

The algorithm starts with setting the initial structural parameters of the BDL and an initial temperature t₀ = 100,000. The temperature drops as $T(i+1)=aT(i)$, where $a=0.85$, $i=0,1,2,3,\dots ,n$, and $n$ is the number of iterations. With the continuous decline of the temperature parameter, the global optimal solution is expected to be obtained^[29]. In terms of the structure parameters, the irradiated diameter of the BDL and the number of concave–convex rings are set as 25 mm and 20, respectively. The values of the radii are defined as variables by setting the minimum radius $t=500\mathrm{\mu m}$. We assume the focal length is 25 mm, and the intensity distribution at the back focal plane is calculated by Eq. (1).

The optimization process is to search for the variable value that satisfies the cost functions, i.e., the FWHM of the main lobe is expected to be no larger than $0.4\lambda$; meanwhile, the sidelobe intensity ratio is less than 20%. The algorithm converges after 88,578 loops; the structural distribution of the THz BDL is illustrated in Fig. 1. The computing time on a laptop with 2.6 GHz i7-10750H, 16 GB RAM, and R2016a MATLAB is about 24 h.

Photosensitive standard resin (PSR) is selected as the 3D printing material. Its transmittance and the refractive index are 82% and 1.66 at 278.6 GHz, respectively^[31]. It is shown in Fig. 1(a) that the yellow concentric rings have a phase modulation of $\pi$, and the rest are 0 by setting the step height $\mathrm{\Delta }h$ at 0.81 mm. In order to facilitate mirror holder clamping, the thickness and edge width of the substrate are set to 1.5 and 2 mm, respectively. The 3D model of the THz BDL is built.

This BDL is fabricated by layer superposition using selective curing printing technique (J850 Pro, Stratasys) with the average thickness of 14 µm of each layer. The liquid tank is filled with liquid PSR. At the beginning of the forming process, the lifting workbench is one layer below the liquid level. The PSR in the scanned area is rapidly cured to a thin sheet when irradiated by a converged laser beam. Once one layer is fabricated, the workbench drops a constant height and cures another layer. The minimum lateral printing size of this printer is 100 µm, which satisfies the accuracy requirements of the BDL working at 278.6 GHz.

The intensity distribution at the back focal plane is simulated by angular spectrum propagation^[25] using MATLAB when the effective area of the BDL is irradiated normally by a 278.6 GHz parallel beam. The pixel pitch is identical with the detector, while the computation window is $2048\mathrm{pixels}\times 2048\mathrm{pixels}$. It is shown in Fig. 2(a) that the central principal maximum is surrounded by several sidelobes. In order to evaluate the diameter of the focal spot, the vertical profile through the beam center is plotted in Fig. 2(b). It is shown that the FWHM of the focal spot is 0.44 mm ($0.4\lambda$), accompanied by two strong sidelobes, indicated by white arrows in Fig. 2(a), which are $\sim 2$ and $\sim 4\mathrm{mm}$ away from the center, respectively. The sidelobe intensity ratio of these two sidelobes is $\sim 18\%$, while ratio values of the other sidelobes are no larger than 10%.

Figure 2(c) shows the intensity distribution recorded by a microbolometer (MICROXCAM-384i-THz, INO) at the back focal plane, with a pixel pitch of $35\mathrm{\mu m}\times 35\mathrm{\mu m}$ featuring $384\mathrm{pixels}\times 288\mathrm{pixels}$ and an exposure time of 50 ms. The THz beam is emitted by an avalanche photodiode (IMMPATT, Terasence) at a wavelength of 1.08 mm (278.6 GHz) with a maximum output power of 26 mW. It is noted that the beam is expanded to 20 mm; therefore, the BDL is not fully irradiated. It is indicated in Fig. 2(c) that there are two detectable sidelobes with intensity ratio of $\sim 22\%$. The FWHM of the actual focal spot, extracted from the curve in Fig. 2(d), is 0.66 mm ($0.6\lambda$), validating that the subwavelength focusing is achieved by the proposed THz BDL. It is also validated that there is a superoscillation effect in this BDL, since the FWHM is smaller than $0.38\lambda /\mathrm{NA}$^[32], which is 0.89 mm (numerical aperture, $\mathrm{NA}=0.462$). When the parallel beam is incident on a BDL with concentric rings, the propagation can be regarded as a coherent superposition of parallel beams with different spatial frequencies. Thereafter, a superoscillation field in the far field is generated by propagating diffracted beams, leading to a subdiffraction focusing spot^[33]. Compared with the simulated beam profile, the inconsistency mainly comes from different irradiated areas of the BDL. Figures 2(e) and 2(f) show simulation results when the irradiated diameter is 20 mm. The FWHM of the main lobe is 0.60 mm ($0.55\lambda$), which is quite close to the experimental value.

The normalized axial intensity distributions of the beam emitted from the BDL are also plotted in Fig. 3. The propagation range is 10 mm to 50 mm, with an axial interval of 0.3 mm. Figure 3(a) shows the simulated distribution when the BDL is fully irradiated like Fig. 2(c). It is shown that the focal length and depth are $\sim 25$ and $\sim 3\mathrm{mm}$, respectively. When the diameter of the incident beam at the BDL is 20 mm, the simulated focal depth shown in Fig. 3(c) matches with the measured value shown in Fig. 3(b), both of which are shortened to $\sim 24\mathrm{mm}$.

The diagram of the experimental setup is illustrated in Fig. 4. To prevent overexposure of the detector, an attenuator with 3% transmittance is placed downstream the source. The divergent beam is collimated by a TCL to form a quasi-parallel beam with a diameter of $\sim 20\mathrm{mm}$. The BDL is placed 49 mm behind the lens and 24 mm in front of the sample. The sample is mounted on a two-axis motorized translation stage (MC600, Zolix). The scanning step is 0.3 mm. The data is recorded by a Golay cell (GC-1D, Tydex) at a distance of 2 mm away from the sample. A chopper is inserted between the TCL and the BDL. The acquisition speed is 0.3 s per pixel.

To quantify the resolution of the THz BDL scanning imaging setup, a 1D resolution test target is fabricated by punching vertical bars with different widths on a stainless-steel plate with a thickness of 2 mm. The scanning region is plotted using a yellow dashed line in Fig. 5(a), in which the smallest width and spacing of the bars are both 0.8 mm, and the largest width and spacing are 2.5 mm. Figure 5(b) illustrates the imaging result using BDL without extra image processing on purpose. The targets with the width of 1, 1.2, 1.5, 1.8, 2, and 2.5 mm can be recognized. Although the edge diffraction from the hollow bars has a negative effect on image quality, there is a clear trend that the signal-to-noise ratio (SNR) increases with width extension of the bars.

For comparison, a TCL (TPX-D25, 1-F25, $\mathrm{NA}=0.447$) is used instead of the BDL at the same location. The FWHM of its focal spot is $\sim 2\mathrm{mm}$. Figure 5(c) presents the corresponding imaging result. The minimum resolution using the TCL is 1.5 mm. Meanwhile, the SNR is lower, and the distortion is more severe than using the BDL. The comparison of the line scan plotted in Fig. 5(d) reveals that using a BDL has higher resolution and quality than using the conventional TCL. The mean error between the measured width using the BDL and the actual width is 15%, while the mean error of the TCL is 25%.

Another sample is a baseplate embossed with the letter “H” with a height of 2 mm. More size parameters are denoted in Fig. 6(a). It is fabricated by 3D printing using visually opaque resin, the transmittance and the refractive index of which are 88% and 1.65 at 278.6 GHz, respectively^[31]. The scanning area is $30\mathrm{mm}\times 30\mathrm{mm}$, as plotted in the mahogany dashed line in Fig. 6(a), with a mechanical translation step of 0.3 mm. A lid made of the same material is 5 mm above the letter, with a thickness of 2 mm.

The imaging results of the sample without or with the lid scanned by the TCL and BDL are shown in Figs. 6(b)–6(e), respectively. In each figure, the blue curve in the subplot corresponds to the position plotted by the white dashed line, indicating the edge undulation. The comparison shows that the letter in Fig. 6(b) has a clearer silhouette than that in Fig. 6(c). In addition, the area of the letter in Fig. 6(b) is $\sim 5\%$ larger than the actual size, and the area of the letter in Fig. 6(c) is $\sim 8\%$ larger than the actual size. The presence of the lid causes a more severe deviation between the size extracted from image and the actual size. When the sample is scanned by the BDL, the area expansion factor is $\sim 2\%$ in Fig. 6(d) and $\sim 3\%$ in Fig. 6(e). By virtue of the smaller focal point and longer focal depth, the BDL also has better imaging quality for the covered object than a conventional TCL made by TPX.

In this Letter, a binary diffractive lens achieving subwavelength focusing is proposed at 278.6 GHz. The diameter of the focal spot is modulated by sets of concentric rings with constant step height and optimized radii; thus, this component can be easily fabricated by 3D printing. The FWHM of the focal spot is 0.65 mm, which varies with the irradiated area of the BDL. A THz scanning imaging setup is constructed based on this low-cost lens, achieving 1 mm ($0.9\lambda$) lateral resolution. The BDL is capable of both converging a focal spot beyond the diffraction limit, which is smaller than that of the TCL, and achieving longer focal length than near-field components, e.g., solid immersion lenses and metalenses. In addition, the fabrication cost of 3D printing is relatively lower than that of the other fabrication approaches, e.g., laser direct-writing^[18], chemical ablation^[19], and photolithography^[20]. By virtue of flexible 3D printing, this phase-type BDL has broad application prospects at the band less than 1 THz, providing a flexible approach for THz beam shaping and phase modulation.

The diameter of this BDL is mainly determined by two factors. First, its efficient area depends on the diameter of the illumination THz beam. To expand a bigger beam, a larger lens upstream from the BDL is required with incremental absorption and transmission loss to the THz radiation. Meanwhile, a reasonable margin of the BDL, 20% in our case, is also very helpful for the alignment in this visually invisible illumination scenario. Second, the computing time based on simulated annealing algorithm would increase with the number of rings. It is feasible to fabricate a bigger BDL, which may further enhance the diffraction efficiency and converge a smaller spot. However, smaller FWHM refers to the main lobe with lower intensity; accompanying sidelobes with stronger intensity would degrade the imaging quality. It is noted that the TCL has higher THz transmission efficiency ($\sim 84\%$) and focusing efficiency ($\sim 71\%$) than the BDL, which are $\sim 72\%$ and $\sim 53\%$, respectively. The decrements are mainly due to lower transmittance of PSR than TPX as well as the existence of sidelobes.

Further study will be conducted to enhance the performance of the BDL and the imaging setup. First, it is feasible to suppress sidelobes by optimizing the BDL structure using high efficiency design approaches, e.g., a deep-learning framework. Second, the nonuniform THz beam is also a main factor degrading the convergence ability of the BDL. It is necessary to develop THz spatial filters based on extraordinary optical transmission to polish the incident beam. Third, for raster scanning when using a cell detector, it is difficult to investigate dynamic phenomena beyond the potential mechanical translation error. A fast scanning system constituted by a 2D galvanometer and array detector is expected to match well with THz BDL scanning.