Main

The past decades witnessed substantial developments in terahertz imaging technologies and their unique utilization in a broad range of applications. The relatively low energy of terahertz photons and their high penetration through many non-conductive materials make terahertz radiation promising for non-destructive biomedical imaging1,2,3; security screening4,5,6; quality control of pharmaceutical7,8,9, industrial10,11,12 and agricultural products13,14,15; and cultural heritage conservation16,17,18.

Despite its great potential, the usage of terahertz imaging in real-world applications has been limited by the slow speed, large size, high cost and complexity of existing imaging systems. Terahertz time-domain spectroscopy (THz-TDS) systems have been instrumental for various imaging applications19,20,21,22,23,24,25,26,27,28,29 since they provide the time-resolved response of objects with a sub-picosecond temporal resolution, which yields both amplitude and phase information over a broad terahertz-frequency range. However, due to the single-pixel nature of existing THz-TDS systems, image data are acquired by raster scanning either the object or the imaging system, leading to very slow imaging speeds and bulky, complex setups.

To avoid raster scanning, electro-optic processes in nonlinear crystals have been employed to convert the object’s terahertz amplitude and phase information to the near-infrared regime and acquire them with an optical camera22,23,24,25. However, due to the nonlinear nature of the wavelength conversion process, these imaging systems generally require bulky and expensive high-energy lasers and provide very low signal-to-noise ratio (SNR) levels (Supplementary Table 1). Another way to avoid raster scanning is to encode the terahertz radiation pattern interacting with the object through a time-varying spatial modulator and reconstruct the image using the spatial pattern information26,27,28,29,30. This approach enabled faster and more robust terahertz imaging systems by eliminating the mechanical stages used for raster scanning. However, these imaging systems require multiple time intervals to encode the terahertz radiation pattern with different spatial distributions, limiting their speed, especially when acquiring image data over a broad terahertz bandwidth. In addition, the requirement for a time-varying spatial modulator adds to the size, complexity and cost of the imaging system.

Further advancements in terahertz imaging technology would substantially benefit from terahertz focal-plane arrays (THz-FPAs) that can directly provide the spatial, ultrafast temporal, spectral, amplitude and phase information of the object simultaneously. Existing terahertz detector arrays based on field-effect transistors31 and microbolometers32 do not provide time-resolved and frequency-resolved image data and lack phase information. It has been shown that phase information can be recovered through digital holography when using these detector arrays in an interferometric setup33,34,35. However, the scope of such interferometric setups is restricted to objects that possess axially uniform refractive indexes and does not apply to more general multilayered structures.

Overall, imaging science in the terahertz part of the electromagnetic spectrum has been lacking FPAs that can directly provide the amplitude and phase information of samples over a large spectral bandwidth and with an ultrafast temporal response. Such capabilities, if made possible in the same FPA, would open up various new applications, for example, difficult-to-see phase-only objects that are weakly scattering could be imaged and sensed in a snapshot, potentially revealing their conformational changes at ultrafast timescales. As another opportunity, broadband spectral features that are acquired from the objects could also reveal their unique spectral signatures, and when this is combined with the amplitude and phase images and the ultrafast temporal response of the samples, it could provide unprecedented channels of information for automated three-dimensional analysis and quantification of samples at high throughput. These powerful information channels, once bundled together, would help us fully exploit all the advantageous features of terahertz waves and substantially improve the resolution and throughput of techniques that are used for imaging, sensing and inspecting materials, objects, scenes and chemical processes in the terahertz-frequency range.

Motivated by these pressing needs, here we present the first THz-FPA that can directly provide the spatial amplitude and phase distributions, along with the ultrafast temporal and spectral information of an imaged object. It consists of a two-dimensional (2D) array of 283,500 plasmonic nanoantennas36,37,38,39 engineered to detect broadband terahertz radiation with a high SNR when used in a THz-TDS system. To simplify data readout from the THz-FPA, these plasmonic nanoantennas are grouped into 7 × 9 clusters, and the collective response of all the nanoantenna clusters is electronically captured at each temporal point to simultaneously resolve their time-domain response. The amplitude and phase responses of the THz-FPA outputs are extracted from the time-domain data over a 3 THz bandwidth. Using an electronic readout, the THz-FPA outputs at each temporal point are captured in 164 μs, exhibiting a 1,000-fold increase in speed compared with the fastest-reported multipixel terahertz time-domain imaging system28, enabling time-domain terahertz video capture at 16 fps (Supplementary Fig. 8 and Supplementary Video 1). Although the use of plasmonic nanoantenna arrays for single-pixel terahertz detection has been previously demonstrated36,37,38,39, terahertz detector arrays based on plasmonic nanoantennas were not realized before, preventing the direct detection of spatial amplitude and phase distributions, as well as the ultrafast temporal and spectral information of an imaged object.

In this work, we utilize the broadband nature of the spatial amplitude and phase data captured by our THz-FPA to realize pixel super-resolution (PSR) using a convolutional neural network (CNN) trained with deep learning. As a first proof of concept, we successfully imaged and super-resolved etched patterns in a high-resistivity silicon substrate and reconstructed both shape and depth of these patterns with an effective number of pixels that exceeds 1 kilopixels, digitally increasing the space–bandwidth product of the THz-FPA by 16 times. This PSR-enhanced THz-FPA would be transformative for building high-throughput and high-resolution imaging and inspection systems that benefit from the unique spectral features of terahertz waves, and would open up a plethora of new applications in biomedical imaging, defense/security screening, agriculture and quality control, among many others.

Results

Figure 1 shows the schematic and operation principles of the terahertz imaging system. It consists of a THz-FPA composed of a 2D array of 283,500 plasmonic nanoantennas fabricated on a low-temperature-grown GaAs (LT-GaAs) substrate, where the plasmonic nanoantennas serve as photoconductive terahertz detector elements. The geometry of the nanoantennas is chosen to obtain a strong spatial overlap between the received terahertz radiation and optical pump beam, despite the large difference between their wavelengths, to achieve high detection sensitivity over a broad terahertz-frequency range (Supplementary Fig. 1). Unlike previously used plasmonic nanoantenna arrays for single-pixel terahertz detection36,37,38,39, which utilized deeply subwavelength nanoantenna lengths, the designed THz-FPA benefits from a wavelength-range nanoantenna length, which provides substantially stronger terahertz electric field enhancement in the device active area. Also, unlike the earlier architectures employed in single-pixel terahertz detection36,37,38,39, the geometry of the plasmonic nanoantennas and shadow electrodes are chosen to accommodate the connection traces of the THz-FPA and maintain a strong overlap between the optical and terahertz beams. To simplify data readout from the THz-FPA, the plasmonic nanoantennas are grouped into 7 × 9 clusters (that is, pixels) and the collective response of all the clusters is captured using a custom-made readout circuit packaged with the FPA chip (Supplementary Fig. 3 and Methods). Using this configuration, ultrafast temporal and spectral information of all the FPA pixels are resolved with more than a 60 dB SNR and a 3 THz detection bandwidth when used in a THz-TDS setup (Supplementary Fig. 4 and Methods). Although this clustering reduces the number of physically accessible pixels of the FPA to 63, the distributed nature of the 283,500 plasmonic nanoantennas positioned uniformly across the entire image plane—when combined with their high SNR and large bandwidth—allowed us to perform PSR using a trained CNN, digitally increasing the space–bandwidth product of the THz-FPA by 16-fold and achieving an effective pixel number of >1 kilopixels. Furthermore, the high sensitivity and broad bandwidth of the THz-FPA enable resolving spectral images of objects over a 2.5 THz bandwidth without using raster scanning or spatial modulation of the terahertz beam (Supplementary Fig. 6).

Fig. 1: Plasmonic photoconductive THz-FPA with PSR.
figure 1

The FPA captures hyperspectral, time-resolved terahertz images of an object and a PSR neural network processes these images to reconstruct a higher-resolution image with a 16-fold larger number of effective pixels. The red and blue beams represent the optical pump and terahertz beams, respectively. λm = 398 μm denotes the median wavelength of the terahertz-frequency band of interest. The inset shows the optical microscopy images (left and middle) of the THz-FPA and a scanning electron microscopy image of the plasmonic nanoantennas (right).

As a proof of concept, we experimentally realized PSR and digitally increased the effective pixel count of our THz-FPA by more than an order of magnitude. For this, we utilized a convolutional deep neural network that was trained on the basis of the experimental image data. Figure 2a illustrates the architecture of the convolutional deep neural network used in this study to process the information recorded by the THz-FPA (Methods). To train the PSR network with our THz-FPA, we imaged spatially structured objects comprising etched patterns in a high-resistivity silicon substrate (Supplementary Fig. 7). Two different experimental setups were designed and tested to image various patterned objects at a distance of 1.1 mm and 20 cm (far field) from the THz-FPA (Supplementary Fig. 5). The same PSR framework was applied to super-resolve the images that are projected onto the THz-FPA in these setups, making our results broadly applicable to different lens-based and lens-free imaging systems.

Fig. 2: PSR-enhanced THz-FPA imaging results.
figure 2

a, The presented imaging system relies on the high SNR and broad detection bandwidth of our THz-FPA and exploits the spatiotemporal information collected in a single shot without any time multiplexing or raster scanning, to enhance the resolution of the imaging system beyond the limit dictated by the physical size of the pixels with the help of a trained CNN. b, Ground-truth thickness images of the objects used in our experiments. Blind testing set A/B contains 54/14 measurements collected based on 27/8 different object patterns. Compared with the training set, the samples in set A differ in terms of their thickness contrasts, whereas the samples in set B are unique with regard to their 2D patterns and thickness contrasts. c, Super-resolved images using the PSR-enhanced THz-FPA. The PSR neural network provides a 16-fold increase in the number of effective pixels, revealing >1 kilopixels in each image. d, Relative spectral phase slope of the objects computed based on the slope of unwrapped spectral phase distribution. e, Amplitude channel at the median frequency of fm = 0.754 THz (λm = 398 μm).

First, we fabricated 300 spatially structured silicon samples, each corresponding to a unique combination of 75 designed patterns and 4 etch depths (Δt = 10, 20, 30 and 40 μm). The objects were imaged at an axial distance of 1.1 mm from the active area of the THz-FPA using a high-refractive-index semiconductor substrate, leading to a Fresnel number of >1. Therefore, one of the main sources of spatial resolution loss for these imaging experiments is the pixelation at the THz-FPA. The resulting measurements from these samples were divided into two sets for the training and blind testing of our PSR-enhanced THz-FPA. The etched depths of all samples, \({\Delta }_{\rm{t}}^{\rm{m}}\), were measured using a profilometer to establish the ground-truth images during the training and testing of our PSR deep neural network.

Our PSR neural network aims to benefit from the large SNR and broad bandwidth of each FPA pixel to achieve super-resolution, and therefore, it is designed to take the raw spectral amplitude and phase components of the THz-FPA output (Fig. 2a). The spectral amplitude components of each FPA pixel at N different frequencies within [fmin, fmax] are used as N input channels to the PSR network. We also enriched this input information by concatenating a thickness (that is, optical path length) image estimate that is computed based on the slope of the unwrapped phase detected within [fmin, fmax]. As a result, we have a total of N + 1 information channels that feed the PSR network, which was trained using error backpropagation with a structural loss term calculated against the ground-truth images. This is a one-time training process, which utilizes the high SNR and multispectral output of our THz-FPA to digitally achieve PSR and increase the space–bandwidth product of our FPA. Although the detection bandwidth of the THz-FPA extends from 0.1 to 3.0 THz, in our training and testing, we took fmin and fmax as 0.5 and 1.0 THz, respectively, avoiding some of the water vapour absorption lines, setting the number of amplitude channels at the input of the PSR network to N = 41. As a result of this choice, the PSR network was trained to process the spatial information contained in the N + 1 = 42 input channels detected by our THz-FPA to perform PSR.

To train and test the CNN shown in Fig. 2a, we used an experimental dataset of 343 raw image measurements conducted on the 300 samples, with some of the samples imaged more than once to demonstrate the repeatability of the imaging system. This image data were partitioned in a way that 68 measurements were preserved for the blind testing of the PSR network, whereas the rest of the measurements were used for training (Supplementary Fig. 15 shows the ground-truth thickness images of the training dataset). Figure 2b depicts the ground-truth images of these 68 test samples that were never used or seen by the PSR network before, and the corresponding super-resolved images are reported in Fig. 2c. Visual comparison of the ground-truth thickness images and the super-resolved PSR network output images reconstructed by taking advantage of the broad bandwidth operation of our multispectral THz-FPA reveals a close match for all the test samples regardless of their material thickness and etch depth. These results also demonstrate resolving small-scale target objects, including patterns with spatial features as small as one-fourth the physical pixel size of our THz-FPA in each lateral direction, which clearly highlights the success of the PSR network. We also performed a quantitative comparison between the ground-truth thickness images shown in Fig. 2b and the neural network’s reconstructions in Fig. 2c using the structural similarity index measure (SSIM) and the peak signal-to-noise ratio (PSNR); our reconstruction results achieved an SSIM of 0.839 ± 0.098 and a PSNR of 16.60 ± 3.67 dB, further revealing the success of our PSR-enhanced THz-FPA.

To characterize the experimentally achieved resolution of our PSR-enhanced THz-FPA, Fig. 3 provides a closer look at the reconstruction performance of the neural network for some grating patterns in our test data with sub-pixel linewidths corresponding to 0.75, 0.50 and 0.25 times the physical pixel size of the FPA. The cross-sections of the ground-truth and PSR network output images with different sub-pixel linewidths (Fig. 3) demonstrate that we can resolve horizontally and vertically oriented thickness variations as small as one-fourth the physical pixel size of the THz-FPA. This indicates a 16-fold increase in the number of effective pixels and a super-resolved total pixel count of 1,008. Despite the deeply subwavelength etch depths of these gratings, the PSR network resolved the 2D structure of these sub-pixel gratings and accurately quantified the thickness variations as shown in the cross-sectional plots (Fig. 3).

Fig. 3: Resolution quantification of the PSR-enhanced THz-FPA.
figure 3

Resolved images of the fabricated resolution test objects are shown on the left. The linewidths and separations are 0.75, 0.50 and 0.25 times the FPA pixel size in both horizontal and vertical directions. The raw images show the depth patterns estimated from the raw data of the THz-FPA. λm = 398 μm denotes the median wavelength of the terahertz-frequency band of interest. The right plots show the averaged depth profiles in the vertical (top three plots) and horizontal (bottom three plots) directions.

To shed more light on the generalization performance of our PSR-enhanced THz-FPA, we trained three additional PSR networks based on different partitioning schemes of the experimental data between the training and testing sets. Supplementary Fig. 9 shows these new tests and their corresponding super-resolved images, all of which successfully demonstrate the generalization capacity of our PSR approach under different training and testing sets.

To further quantify the imaging performance of our PSR-enhanced THz-FPA, in Fig. 4a, we report the depth estimation \({\hat{\varDelta }_{{\rm{t}}}}\) for each of the reconstructed test patterns shown in Fig. 2b, compared with the ground-truth values \({\Delta }_{\rm{t}}^{\rm{m}}\) measured with a profilometer. As desired, \({\hat{\Delta }_{\rm t}}\) inferred by the PSR network output images generally follow the \({\Delta }_{\rm t}^{\rm m}={\hat{\Delta }_{\rm t}}\) line (Fig. 4a, dashed line). Figure 4b–d illustrates a similar analysis for the three additional test datasets (Supplementary Fig. 9), where the inference results also follow the \({\Delta }_{\rm t}^{\rm m}={\hat{\Delta }_{\rm t}}\) line (Fig. 4b–d, dashed lines) without any major errors or outliers, demonstrating the success of the presented PSR-enhanced THz-FPA in achieving both super-resolution and material thickness estimation, independent of the training and testing datasets. This conclusion is also supported by the SSIM and PSNR values shown in Fig. 4e,g. We also analysed the SSIM and PSNR of the reconstructed images across all the test sets using the percentage of thickness estimation error, defined as \(\frac{\left|{\Delta }_{\rm t}^{\rm m}-{\hat{\Delta }_{\rm t}}\right|}{{\Delta }_{\rm t}^{\rm m}}\times 100\) (Fig. 4f). Based on this analysis, the mean error in \({\hat{\Delta }_{\rm t}}\) inferred by our PSR network is found to be 10.73%; for example, for the test samples fabricated with Δt = 10 µm, this mean error corresponds to ~1 μm, which is deeply subwavelength.

Fig. 4: Quantification of the image reconstruction quality and thickness contrast prediction accuracy.
figure 4

a, \({\hat{\Delta }_{\rm t}}\) versus \({\Delta }_{\rm t}^{\rm m}\) for the test set shown in Fig. 2. The dashed black line represents \({\hat{\Delta }_{\rm t}}=\) \({\Delta }_{\rm t}^{\rm m}\). bd, Same as a, except that the reported \({\hat{\Delta }_{\rm t}}\) and \({\Delta }_{\rm t}^{\rm m}\) represent the estimated and true (ground-truth) values of material thickness contrasts for the test image sets, namely, set 2A and set B (b), set 3A and set B (c) and set 4A and set B (d) (Supplementary Fig. 5a–c). e, Violin plots of SSIM (left) and PSNR (right) values representing the image reconstruction quality achieved by the PSR networks. The horizontal line in each violin plot shows the mean value. The box plot inside each violin plot indicates the interquartile range, and the dot denotes the median value. Nm represents the sample size for each test set. The only difference between these PSR networks is their partitioning of the training and testing sets (Methods). f, Distribution of the SSIM (left) and PSNR (right) of all the object patterns in the test sets, namely, set A, set B, set 2A, set 3A and set 4A as a function of the percentage of thickness estimation error, that is, \(\frac{\left|{\hat{\Delta }_{\rm t}}{-\Delta }_{\rm t}^{\rm m}\right|}{{\Delta }_{\rm t}^{\rm m}}\times 100\). g, Violin plots of SSIM (left) and PSNR (right) values along with the corresponding standard deviations representing the image reconstruction quality achieved by the PSR networks for the object patterns included in the test sets of set A, set B, set 2A, set 3A and set 4A as a function of Δt. Nm denotes the sample size for each test set.

To further assess the generalization performance of our PSR-enhanced THz-FPA, we designed, fabricated and imaged various objects with more complex and irregular spatial structures than the ones used for the training, including hand-written numbers and letters, and verified the robustness of our PSR algorithm in resolving both depth and shape of these objects. These newly fabricated objects were blindly tested with the same PSR neural network used in Fig. 2, and no patterns with similar shapes were seen by the neural network during its training. The depth and shape of the patterns were successfully reconstructed with an average PSNR of 12 dB and SSIM of 0.64, confirming the generalization performance of the presented PSR-enhanced THz-FPA (Supplementary Fig. 11). In addition, we numerically tested the same PSR-enhanced THz-FPA framework for resolving images of objects with non-binary depths and weaker contrast. Supplementary Fig. 16 shows these new objects with variable depths/thicknesses and their corresponding super-resolved images; our results reveal the successful reconstruction of these test object images with an average PSNR of 29.4 dB and SSIM of 0.93.

To demonstrate the capability of our PSR-enhanced THz-FPA for imaging far-field objects, we imaged several spatially structured silicon samples with four different etch depths (Δt= 10, 20, 30 and 40 μm), placed at an axial distance of 20 cm from the THz-FPA surface (Supplementary Fig. 5). We constructed an experimental image dataset containing 151 far-field measurements, selected 16 measurements corresponding to 16 unique samples for blind testing and used the remaining 135 measurements to train the PSR network (Supplementary Fig. 20 shows the ground-truth thickness images of the training dataset). In our training and testing, we took fmin and fmax as 1.3 THz and 2.1 THz, respectively, avoiding some of the water vapour absorption lines, setting the number of amplitude channels at the input of the PSR network to N = 52. The ground-truth images of the 16 test samples representing our blind testing set and the corresponding super-resolved images are shown in Fig. 5a,b, respectively. Visual comparison of the ground-truth thickness images and the super-resolved images indicates a close match for all the test samples regardless of their thickness and etch depth. Our reconstruction results achieved an average SSIM of 0.78 and a PSNR of 14.5 dB, demonstrating the success of our PSR-enhanced THz-FPA imaging system in far-field measurements.

Fig. 5: Far-field PSR-enhanced THz-FPA imaging results.
figure 5

a, Ground-truth thickness images of the objects used in our experiments. b, Super-resolved images using the PSR-enhanced THz-FPA. c, Relative spectral phase slope of the objects computed based on the slope of unwrapped spectral phase distribution. d, Amplitude channel at the median frequency of fm = 1.6 THz.

Finally, to experimentally assess the resilience of our PSR-enhanced THz-FPA imaging system to potential object misalignments, we imaged test objects (1) at different axial distances from the THz-FPA and (2) with different rotation angles within the sample field of view, and demonstrated the successful reconstruction of both shape and depth of these purposely misaligned objects despite the fact that the network training was free from such large degrees of misalignments (Supplementary Figs. 14, 18 and 19).

Discussion

In this work, we present one of the first THz-FPAs comprising 283,500 plasmonic photoconductive nanoantennas, capable of simultaneously providing the spatial amplitude and phase distributions as well as the ultrafast temporal and spectral information of an imaged object. Our PSR-enhanced THz-FPA exploits its multispectral operation to increase the spatial resolution and the effective number of pixels in the final reconstructions compared with the physical pixel size of the FPA. The presented deep learning-driven terahertz imaging system uses the amplitude and phase information of the detected spectral components of an imaged object to enhance the spatial resolution by 4-fold in each lateral direction, increasing the space–bandwidth product of the FPA by 16-fold (4 × 4), with >1 kilopixels in each reconstructed image. In addition, our PSR-enhanced THz-FPA has been shown to quantify deeply subwavelength material thickness variations as small as 10 µm, with a mean thickness accuracy of ~1 µm.

To further emphasize the importance of the multispectral operation of the THz-FPA, we compared the performance of our PSR framework with multispectral input data against the use of single-frequency input data. For this purpose, we trained six image reconstruction deep neural networks that use the amplitude and phase outputs of the THz-FPA at a single frequency of 0.50, 0.75, 1.00, 1.25, 1.60 and 2.00 THz. Supplementary Fig. 10 shows the reconstructed images of each one of these single-frequency PSR networks along with their corresponding SSIM and PSNR values. Although these single-frequency systems can reveal the images of the input objects to some extent, the reconstructed images suffer from severe artifacts and distortions compared with our multispectral PSR results (Fig. 2). These results testify that a single-frequency or narrowband THz-FPA capable of providing both amplitude and phase information cannot match the super-resolution capabilities offered by our multispectral THz-FPA. It should also be noted that the existing detector arrays based on field-effect transistors31 and microbolometers32 do not provide frequency-resolved image data and cannot obtain phase information, which would prevent their performance enhancement with the discussed PSR framework.

Our neural network-based PSR results demonstrated strong generalization to new imaging conditions (for example, object rotations and axial defocusing) and new types of sample never used in the training before (Supplementary Figs. 11, 14, 16, 18 and 19). The quality of the resolved images using the reported PSR framework could be further improved with additional training data and unique objects, which would also allow us to make use of deeper neural network architectures without overfitting. Additionally, we did not use a generative model trained through an adversarial game between a generator and a discriminator network; instead, we solely relied on a structural loss function in the form of mean squared error. With sufficiently large training data and carefully balanced loss terms using a linear mix between a structural and adversarial loss, the image reconstruction quality reported in this study could be further improved using, for example, generative adversarial networks40,41.

The presented multispectral super-resolution framework used for imaging phase objects demonstrates one potential application of our plasmonic photoconductive THz-FPA as its proof of concept. To demonstrate that the imaging capabilities of this THz-FPA are not limited to super-resolving phase-only objects through a neural-network-based deep learning algorithm, we also captured and super-resolved a video (at 16 fps) of water flow in three adjacent plastic pipes (each with an inner diameter of 250 μm) by using a completely different algorithm based on holography (Supplementary Fig. 21 and Supplementary Video 1). Since water strongly absorbs terahertz waves, the channels that carry water form an opaque object that is dynamically evolving through the flow, in contrast to the transmissive phase objects that were used earlier in our experiments.

Based on all the proof-of-concept experimental results reported in this work, we believe that the unique capabilities of our plasmonic photoconductive THz-FPA enable the high-speed acquisition of spatial, ultrafast temporal, spectral, amplitude and phase information of the imaged object; this diverse set of information acquired through our THz-FPA could be used in various applications by, for example, utilizing the time-of-flight terahertz pulses to reconstruct three-dimensional images of multilayered objects42, using the spectroscopic signatures to identify chemicals19,20,21 and integrating it with diffractive optical networks for feature detection and object classification43,44. For example, the demonstrated THz-FPA offers a powerful and much needed solution for a long-lasting problem in battery manufacturing: non-destructive quality assessment of battery electrodes during manufacturing. By extracting volumetric structural information of battery electrodes, which are layered structures of optically opaque coatings on metallic plates, our THz-FPA enables the real-time and high-throughput detection of various types of hidden defect such as divots, agglomerates and metal contamination, as well as non-uniformities in coating thickness (Fig. 6). These defects and non-uniformities drastically change the electrochemical performance of batteries, diminishing their performance by aggravating cycle efficiency, lowering discharge capacity, shortening the life span and causing safety hazards like fire and explosion. Figure 6 illustrates different defective lithium-ion battery electrodes quantitatively imaged by our THz-FPA and the resolved intensity and thickness profiles, revealing structural information of these defective electrodes. To date, no other imaging method has been capable of detecting and quantifying such structural defect information in real time, to the best of our knowledge, which is urgently needed in roll-to-roll battery manufacturing plants considering all the applications and products that heavily rely on high-quality batteries.

Fig. 6: Structural information of defective lithium-ion battery electrodes captured with our THz-FPA.
figure 6

a, Optical camera image of a lithium-ion battery electrode (left), showing a slightly darker region around an agglomerate under a strong flashlight. The top view of the terahertz intensity (middle) and coating thickness (right) profiles captured by the THz-FPA clearly indicate the presence of the agglomerate, which causes terahertz scattering and thickness variations around the defective region. A uniform coating thickness of 100 μm is measured for the healthy regions of the battery electrode. b, Optical camera images of different lithium-ion battery electrodes (left), showing small changes in the darkness of the defective regions (agglomerates and divots) under a strong flashlight. The top view of the terahertz intensity (middle) and cross-sectional (right) profiles captured by the THz-FPA clearly indicate the presence of defects, which causes terahertz scattering and variations in coating thickness around the defective regions. A uniform coating thickness of 100 μm is measured for the healthy regions of the battery electrodes.

The demonstrated THz-FPA can be further advanced to provide a larger number of pixels, larger field of view, higher SNR, higher bandwidth and faster image acquisition. The use of photoconductive nanoantennas integrated with a plasmonic cavity offers considerably higher SNR and bandwidth at a much lower optical power level45, and could enable larger-pixel-count FPAs and larger-field-of-view imaging systems without any degradation in the image quality. The speed of the demonstrated THz-FPA is still limited by the utilized sequential readout electronic system. Integration of the FPA with a 2D parallel readout integrated circuits would substantially increase the imaging speed, SNR and bandwidth. These advancements could bring us much closer to realizing the untapped potential of the terahertz spectrum for numerous applications.

Methods

THz-FPA fabrication process

The THz-FPA was fabricated on an LT-GaAs substrate with a carrier lifetime of 0.3 ps (Fig. 1, inset). Supplementary Fig. 2 shows the fabrication process flow. The plasmonic nanoantennas were patterned using electron-beam lithography, followed by the deposition of 3/47 nm Ti/Au and lift-off. The dipole-shaped nanoantennas have a 100 nm width and 180 nm periodicity along the z axis and have an arm length and tip-to-tip gap size of 20 μm and 500 nm, respectively, along the x axis. Output traces and contact pads were patterned by photolithography, followed by 20/180 nm Ti/Au deposition and lift-off. A 290-nm-thick silicon nitride anti-reflection coating was deposited by plasma-enhanced chemical vapour deposition. Shadow metals were formed on the anti-reflection coating by photolithography, 10/90 nm Ti/Au deposition and lift-off. The fabricated THz-FPA consists of 283,500 plasmonic nanoantennas covering a 2.4 mm × 1.7 mm area.

Terahertz time-domain imaging setup

A Ti:sapphire laser (Coherent Mira 900) was used to generate optical pulses with an 800 nm central wavelength, 135 fs pulse width and 76 MHz repetition rate. The optical beam was split into two branches. One branch pumped a plasmonic photoconductive terahertz emitter46 to generate terahertz pulses. The other branch pumped the THz-FPA. The optical pump power incident on the terahertz emitter and FPA was 660 and 500 mW, respectively. For the experimental setup used for Fig. 2b–e, two parabolic mirrors were used to collimate and focus the generated terahertz radiation onto the active area of the THz-FPA after interacting with the imaged object. The objects were axially placed ~1.1 mm away from the active area of the THz-FPA using a high-refractive-index semiconductor substrate (Supplementary Fig. 5a). For the experimental setup shown in Fig. 5, the imaged objects were placed at a 20 cm axial distance from the active area of the THz-FPA. One parabolic mirror and one terahertz objective lens (TeraLens, Lytid) were used to collimate and focus the terahertz radiation onto the THz-FPA after interacting with the imaged objects and providing a demagnification factor of ~2.76 (Supplementary Fig. 5b). An optical delay stage was used to vary the optical path difference between the optical pump and terahertz beams incident on the THz-FPA. The time-domain terahertz electric field was obtained by measuring the photocurrent of the FPA pixels as a function of the optical delay. The frequency-dependent amplitude and phase of the terahertz signal were calculated by taking the Fourier transform of the time-domain data.

THz-FPA data acquisition system

A custom-made readout circuit was built for data readout from the THz-FPA. The circuit consists of a field-programmable gate array development board (Basys 3, Digilent) controlling four 16-channel multiplexers (ADG1206, Analog Devices), which sequentially direct the FPA outputs to a transimpedance amplifier (DLPCA-200, FEMTO). The transimpedance amplifier converts the FPA photocurrent outputs into voltage signals, which are then sampled by a lock-in amplifier (MFLI, Zurich Instruments) and routed to a computer for post-processing. The THz-FPA outputs at each temporal point are captured in 164 μs.

Phase object fabrication process

The phase objects were fabricated on a high-resistivity silicon substrate. Supplementary Fig. 7 shows the fabrication process flow. A thin layer of silicon dioxide was first grown on the silicon substrate using thermal oxidation. Then, the object patterns were transferred onto the silicon dioxide layer by photolithography and silicon dioxide etching. Deep reactive ion etching was used to create the trenches in the silicon substrate, with the silicon dioxide layer being the etch mask. The depth of the silicon trenches was controlled by the etching time. Finally, the silicon dioxide layer was removed through buffered oxide etching. The depth of the silicon trenches was measured by a Dektak 6M profilometer, which provides a 1 nm depth resolution.

Data preparation for deep neural network training and testing

The generation of ground-truth images was automated using a MATLAB R2020a code that takes the etch mask (.dxf file) of the fabricated objects and outputs their thickness map with a resolution of 10 μm in both x and z directions, based on morphological image processing operations47. Since the relative position of the objects with respect to the active area of the THz-FPA is subject to change during the experiments, we implemented an algorithm that estimates the 2D relative location of the objects and accordingly shifts the high-resolution ground-truth images. For the estimation of the object location, we followed a two-step procedure. In the first step, the high-resolution ground-truth images were downsampled to a Cartesian grid defined by the physical pixel size of our FPA, that is, 240 μm × 270 μm. For each spectral component detected within the 0.5–1.0 THz range, we computed the normalized cross-correlation function between the detected amplitude channels and the downsampled, low-resolution ground-truth model. The spectral component exhibiting the highest correlation coefficient was selected as the reference for the second fine adjustment step, and the location of the peak correlation coefficient was used for the coarse alignment.

Since image alignment accuracy beyond the physical pixel size was required to accurately train the PSR networks, we adapted a greedy strategy for fine alignment. Each ground-truth high-resolution image was shifted to 24 × 27 different positions with 10 μm steps around its coarsely estimated location computed in the first step. At each position, the ground truth was downsampled using bilinear interpolation to the physical resolution of the FPA and compared with the amplitude channel of the reference spectral component determined in the first step. Out of 24 × 27 = 648 normalized cross-correlations, we determined the best alignment position as the one providing the highest correlation coefficient, providing us an estimate of the relative object location with respect to the FPA at a spatial sampling rate much higher than the size of the physical pixels. All of these steps constitute a one-time effort used for training our PSR networks.

To super-resolve the images of the spatially structured silicon samples placed at an axial distance of 1.1 mm from the THz-FPA, we constructed a dataset containing 343 measurements; out of these 343 measurements, we reserved 68 measurements corresponding to 35 unique samples for blind testing (Fig. 2b), whereas the remaining 275 measurements were used to train the PSR network. In these 68 measurements constituting the blind testing dataset shown in Fig. 2, there were two major groups of objects labelled as set A with 54 measurements and set B with 14 measurements. The 54 measurements contained in set A were distributed over 27 unique object patterns, and for each of these object patterns, the training set contained measurements produced by the same 2D pattern although fabricated with different etch thickness levels. The samples in set B, on the other hand, did not only differ from the training samples in terms of the etch thickness contrast but their 2D patterns were also unique, not included in the training set. Preserving this training/testing partition ratio (275/68), we also created three additional datasets (Supplementary Fig. 9). Although the samples in set B of these three new blind testing datasets were the same as the one shown in Fig. 2, the samples in set A were different and complementary to each other in terms of material thickness contrasts. To super-resolve the images of the spatially structured silicon samples placed at an axial distance of 20 cm from the plane of the THz-FPA, we constructed an experimental dataset containing 151 far-field measurements, reserved 16 measurements corresponding to 16 unique samples for blind testing and used the remaining 135 measurements to train the PSR neural network.

PSR network architecture

Figure 2a illustrates the architecture of the deep neural network for PSR. In this study, we used a CNN inspired by the upsampling arm of U-Net48. Although we investigated the performance of other CNN architectures proven to provide high-quality super-resolved images, for example, SRResNet49, SRCNN50 and ESPCN51, the architecture shown in Fig. 2a outperformed these architectures based on our THz-FPA results.

Each grey arrow in Fig. 2a represents a convolutional block that consists of three convolutional layers. If we denote the number of channels at the input (output) of a convolutional block c by \({N}_{c}^{\;{\rm{i}}}\) (\({N}_{c}^{\;{\rm{o}}}\)), then the number of channels is decreased at the output by half, that is, \({N}_{c}^{\;{\rm{o}}}=\frac{{N}_{c}^{\;{\rm{i}}}}{2}\). This reduction in the number of channels is achieved by the three convolutional layers in cth block, where the number of channels at the input and output of each of these convolution operations is given as \(\left({N}_{c}^{\;{\rm{i}}},\frac{{N}_{c}^{\;{\rm{i}}}+{N}_{c}^{\;{\rm{o}}}}{2}\right)\), \(\left(\frac{{N}_{c}^{\;{\rm{i}}}+{N}_{c}^{\;{\rm{o}}}}{2},\frac{{N}_{c}^{\;{\rm{i}}}+{N}_{c}^{\;{\rm{o}}}}{2}\right)\) and \(\left(\frac{{N}_{c}^{\;{\rm{i}}}+{N}_{c}^{\;{\rm{o}}}}{2},{N}_{c}^{\;{\rm{o}}}\right)\). In addition, between each convolutional block, there is a 2 × 2 bilinear upsampling operation. We used four convolutional blocks interconnected through three linear upsampling operations. Including the first convolutional layer that takes Nλ-amplitude channels concatenated with the 7 × 9 low-resolution thickness estimate \({{\bf{I}}}_{{\rm{t}}}^{{\rm{LR}}}\) and the output layer (Fig. 2a, blue arrow), our neural network consists of a total of 14 convolutional layers.

To overcome overfitting, we introduced data augmentation as part of the training forward model, in the form of random vertical and horizontal flipping of the input amplitude and thickness estimate channels. In addition, we employed random 90° rotations, with an additional resampling procedure (using bilinear interpolation) due to the non-symmetric structure of our FPA considering the number of pixels in the x and z directions. Although the random flipping and 90° rotations prevented overfitting and improved the generalization performance of the PSR network, to further enhance our super-resolution imaging performance, we adapted a simple yet effective data augmentation technique, namely, MixUp52,53. The MixUp technique creates linear combinations of input and output pairs based on a random variable in the range [0, 1] drawn from a Beta(α, α) distribution. In this study, α was taken as 0.2. With Xq and Xp denoting the input data sensors for the qth and pth experimental measurements, respectively, the MixUp technique created a new input, namely Xqp:

$${X}_{{qp}}=\delta {X}_{q}+(1-\delta ){X}_{p}.$$
(1)

Similarly, the corresponding ground-truth image \({I}_{{qp}}^{{{\rm{SR}}}}\) was generated/augmented:

$${I}_{{qp}}^{{SR}}=\delta {I}_{q}^{{SR}}+(1-\delta ){I}_{p}^{{SR}},$$
(2)

where \({I}_{q}^{{SR}}\) and \({I}_{p}^{{SR}}\) represent the ground-truth images corresponding to the data tensors Xq and Xp, respectively. All the ground-truth images ISR are measured by the Dektak 6M profilometer.

We trained our PSR networks using a structural loss function L based on the mean squared error as

$$L=\frac{1}{{MN}}\mathop{\sum }\limits_{m=1}^{M}\mathop{\sum }\limits_{n=1}^{N}{\left|{I}^{{SR}}\left[m,n\right]-{I}^{{{SR}}^{{\prime} }}\left[m,n\right]\right|}^{2}.$$
(3)

Here \({I}^{{{SR}}^{{\prime} }}\) defines the inference of the PSR network with M × N output pixels.

The image quality metric, PSNR, which was directly related to the mean squared error, was computed as

$${\rm{PSNR}}=20{\log }_{10}\left(\frac{T}{\sqrt{L}}\right),$$
(4)

where T denotes the etch depth of the ground-truth object measured by the Dektak 6M profilometer.

To compute the other image quality metric, SSIM, we used the following equation:

$${\rm{SSIM}}=\frac{1}{W}\mathop{\sum }\limits_{w=1}^{W}\frac{\left(2{\mu }_{{I}_{w}^{\;{\rm{SR}}}}{\mu }_{{{I}_{\;w}^{\;{\rm{SR}}}}^{{\prime} }}+{C}_{1}\right)\left(2{\sigma }_{{I}_{w}^{\;{\rm{SR}}}{{I}_{w}^{\;{\rm{SR}}}}^{{\prime} }}+{C}_{2}\right)}{\left({{\mu }_{{I}_{w}^{\;{\rm{SR}}}}}^{2}+{{\mu }_{{{I}_{w}^{\;{\rm{SR}}}}^{{\prime} }}}^{2}+{C}_{1}\right)\left({{\sigma }_{{I}_{w}^{\;{\rm{SR}}}}}^{2}+{{\sigma }_{{{I}_{w}^{\;{\rm{SR}}}}^{{\prime} }}}^{2}+{C}_{2}\right)},$$
(5)

where the terms \({\mu }_{{I}_{w}^{{SR}}}\) and \({\mu }_{{{I}_{w}^{{SR}}}^{{\prime} }}\) denote the local mean computed for each 8 × 8 window w, scanned over the ground-truth thickness image ISR and the super-resolved image at the output of the PSR network ISR′, respectively. Similarly, \({\sigma }_{{I}_{w}^{{SR}}}\) and \({\sigma }_{{{I}_{w}^{{SR}}}^{{\prime} }}\) denote the corresponding local standard deviations with \({\sigma }_{{I}_{w}^{{SR}}{{I}_{w}^{{SR}}}^{{\prime} }}\) representing the local cross-covariance between the corresponding windows from ISR and ISR′. As suggested in its original implementation54, we convolved both ISR and ISR′ with an 11 × 11 truncated, circularly symmetric Gaussian filter with a standard deviation of 1.5 before computing equation (5). The constant factors in equation (5), that is, C1 and C2, were computed as (K1T)2 and (K2T)2, where T denotes the ground-truth material thickness contrast. The multiplicative terms K1 and K2 were set to be 0.01 and 0.03, respectively, following the original implementation54.