Terahertz (THz) imaging systems are extensively used to determine the chemical and material composition of samples, marking a pivotal driving motivation in the growth of THz science as a distinct discipline [1–5]. However, the limited availability of a wide field-of-view (FOV) and high-resolution imaging system poses a prominent technological constraint in this field of research. Specifically, at 0.3 mm wavelength (equivalent to 1 THz frequency), the far-field spatial resolution of an image is conventionally constrained to the Rayleigh criterion that corresponds to ∼180μm in vacuum, which often results in a limited space-bandwidth product [6] of the currently available THz imaging systems. To address these challenges, THz digital holography in both off-axis [7–11] and in-line [12–16] configurations has been merged with advanced computational processing algorithms such as sub-pixel sampling [14], extrapolation [17], autofocusing [15], iterative phase retrieval [14,15], and compressive sensing [16] and demonstrated improvement in the image resolution and reconstruction quality. However, prior work on THz holography has employed a monochromatic THz source, and the approach requires the incorporation of broadband THz imaging capabilities. Along the same direction, wavefront shaping techniques have recently shown promise in the THz spectral range and demonstrated single-pixel ghost imaging for complex structural objects [18–26]. However, wavefront shaping approaches involve a rather bulky and complex imaging setup due to the limited capabilities of THz wavefront modulation [27–29]. In addition, THz single-pixel ghost imaging modalities generally depend on near-field detection of scattered fields from the object, which restricts their ability to images in far-field conditions. Similarly, ptychography imaging techniques have been introduced to accelerate high-resolution imaging without significant hardware modifications to THz imaging systems [30–32]. Quite interestingly, THz near-field imaging modalities [33–39], including tip-enhanced THz probes [37,40–42], micro-aperture THz probes [43], and air-plasma dynamic aperture [44,45], have been introduced to achieve deep subwavelength resolution by capturing the object-modulated evanescent waves before the light diffraction occurs. While near-field imaging techniques have achieved deep subwavelength spatial resolution, these conventional methods are implemented through a sophisticated experimental approach and still necessitate mechanical raster scanning of an object surface with a relatively low signal-to-noise ratio (SNR), requiring a relatively long image acquisition time.

Conceptually, the resolution of an image is defined by the quantity of high-spatial-frequency components present in the scattered waveform from the sample [46,47]. Moreover, during free-space propagation, high-spatial-frequency components typically diffract at greater angles than low-spatial-frequency components, causing the diffraction pattern to extend beyond the FOV of the detector. Consequently, a finite numerical aperture (NA) of the optical imaging system restricts the resolution by failing to capture high-spatial-frequency components. In such a direction, the synthetic aperture imaging technique has been introduced to overcome the constraints imposed by the limited NA in imaging systems, which relies on capturing multiple diffraction patterns, each containing a limited spatial frequency spectrum, and coherently combining them to exploit a larger spatial frequency passband for complex samples. Considering this, synthetic-aperture-based imaging was first developed at radio wavelength to improve the resolution of a single radio telescope [48,49] and subsequently adopted at optical [50–52], electron microscopy [53], sonar [54], and THz frequencies [13,55].

A significant challenge in the optical domain is the measurement of complex wavefronts, as image sensors are sensitive to the average intensity distribution of the signal. Such limited availability of field-sensitive detectors complicates the direct assessment of phase information or requires phase-retrieval algorithms [56]. In contrast, THz time-domain spectroscopy (TDS) is an established technique, able to fully resolve electric-field oscillations in very broadband pulses, thereby providing full access to the complex spectral field after the Fourier transform. A meaningful scientific question is whether the THz-TDS can be exploited to develop synthetic aperture imaging approaches where the temporal field information can directly contribute to the high-resolution image reconstruction.

In this work, we explore this unique potential by devising a new approach for time-resolved, coherent synthetic aperture imaging. Our methodology, which outperforms the resolution limit in conventional THz imaging, involves synthesizing an expanded spatial frequency bandpass filter scanned across the Fourier space by illuminating the sample with a broadband THz pulse at various angles. In our imaging analysis, we then incorporate a propagation of the temporal traces obtained through TDS detection, facilitating the reconstruction of the full-field response of the sample, including phase and temporal delay contributions. This approach permits us to acquire extensive spatiotemporal information on the light-matter interaction between the THz beam and the sample. After computational reconstruction, our method offers hyperspectral imaging capabilities that enable the characterization of the complex, subwavelength resolution morphology of a sample across the broadband THz spectral range.

We present the proposed THz Fourier synthetic aperture technique implemented using a TDS imaging system. Figure 1(a) outlines our methodology, where the target object is probed by the series of THz pulses incident at several angles. The time-resolved scattered fields from the object are recorded by raster scanning a single-pixel TDS sensor on the image plane. Our imaging framework can be described as follows. A thin semi-transparent object (neglecting the absorption effect), characterized by a coherent transfer function, T˜(x,y;ω), that we seek to recover, is illuminated by a plane wave defined as exp(i(kxax+kyay))f(ω), with the incident wave vectors k→a≡(kxa,kya) and f(ω) the THz broadband pulse spectrum, typically extending over a few THz. The transmitted field in its spatial frequency spectrum at a particular angular frequency ω can be expressed as E˜+(kx−kxa,ky−kya;ω)=Fx→k2D{T˜(x,y;ω)·ei(kxax+kyay)f(ω)},(1)where Fx→k2D{−} represents the 2D spatial Fourier transform and ‘·’ denotes point-wise multiplication. The transmitted field spectrum, E˜+, then propagates through the pupil aperture, which modulates the spatial frequency spectrum by imposing a low-pass filter. Consequently, the final spatiotemporal field at the imaging plane can be written as Eout(x,y;t)=Fω→t−1D⌊Fk→x−2D{E˜+(kx−kxa,ky−kya;ω)·Op(ω)(kx,ky)}⌋,(2)where Op(ω) denotes the pupil aperture due to the pinhole, Fk→x−2D{−} represents 2D spatial inverse Fourier transform, and Fω→t−1D⌊−⌋ denotes the inverse time-Fourier transform. The overall transmitted spatiotemporal field [in Eq. (2)] at the imaging plane can be experimentally obtained using a combination of a photoconductive-antenna-based THz emitter and detector, readily available commercial products. Notably, the captured TDS images could immediately provide low-resolution broadband field images by performing a time-Fourier transform for each spatial point of the recorded spatiotemporal fields.

For a conventional THz imaging system operating at ω and having a finite NA, the pupil aperture is represented as a circular area formed by kx2+ky2<(NAω2πc)2 in the spatial-frequency space. This means that TDS images formed at the image plane contain spatial-frequency content limited by the pupil aperture. Synthetic aperture imaging exploits the fact that the (−kxa,−kya) wavevector offset in the transmitted field images, as shown in Eqs. (1) and (2), directly correlates with the illumination angle. Therefore, tweaking the illumination angle is equivalent to shifting the center of the pupil aperture, thereby enabling the capture of higher-spatial-frequency components that would otherwise be inaccessible within a conventional imaging framework.

For our THz-TDS proposed implementation, Fig. 1(b) illustrates how changing the illumination angles of the TDS pulses affects the shift in the spectrum of the transmitted field at 0.5, 1, and 1.5 THz. A spatial frequency low-pass filter (shaded circular region in the k-space of the object), typically controlled by a pinhole in the experimental scenarios, further constrains the measurable portion of the spectrum, resulting in a low-resolution and distorted image. A notable difference with optical implementation is that due to their very large spectral bandwidth, the effective size of the pupil aperture is impacted by the frequencies present in the THz pulse, which determines the variable aperture throughout the pulse spectrum (see Visualization 1). Figure 1(c) shows the measured field amplitude images corresponding to the spatial spectrum encompassing the pupil aperture, as shown in Fig. 1(b). It is worth noting that the hyperspectral field images, shown in Fig. 1(c), are directly acquired from the time-resolved measurements as our approach incorporates broadband THz pulse interaction with the sample.

As a next step, a set of spatiotemporal fields is recorded by varying illumination angles of the input THz pulses, and spatial frequency spectra of these transmitted fields are coherently stitched together. This implies that, by scanning the confined pupil aperture in Fourier space, following the spiral pattern shown in Fig. 1(b), we could synthesize a larger spatial frequency passband for inhomogeneous samples. To reconstruct the complex transfer function of the imaging object, we propose to implement a least-squares minimization problem that relies on reducing the error between the measured field E˜outa(x,y;ω) and the expected field ξ˜a predicted by the forward model outlined in Eqs. (1) and (2). In this scenario, we define the fitness function at ω as follows: F(T˜r)=argminT˜r∑aA∑x∑y|E˜outa(x,y;ω)−ξ˜a{T˜r(x,y;ω)}|2,(3)where a corresponds to the input beam having spatial wave vector k→a≡(kxa,kya), A is the total number of illumination angles, and ξ˜a{−} refers to the forward model operator that predicts the spatial field distribution on the image plane after illuminating the sample with the THz pulse of input wavevector k→a. The convex optimization problem defined in Eq. (3) is further solved by employing a Nesterov accelerated gradient-descent algorithm [57], which iteratively recovers the transfer function of the imaging object. We emphasize that the fitness function described in Eq. (3) is field-sensitive, and when minimized, it can simultaneously optimize the amplitude and absolute phase of multispectral images. This gives advantages over time-averaged and intensity-based detection where access to the absolute phase with arbitrary precision is highly challenging, and waveform synthesis requires prior knowledge of the absolute phase profile of the input waveform (e.g., in the case of a transform-limited pulse). Finally, to evaluate the reliability of our numerical approach in the presence of experimental noise, we conducted the numerical analysis under the assumption that the THz source is providing 1 nJ transform-limited pulse of duration 250 fs with 30 dB SNR measured at the detection TDS sensor. Based on this assumption, we added a white-noise component compatible with the experimental scenarios. Such relatively low noise can be readily attained in experiments with many commercial time-domain systems. Moreover, in our simulation, a 100mm2 FOV is raster scanned at 50 μm resolution by the single-pixel TDS detector, mimicking experimental implementations of the THz imaging system [58]. The proposed approach is experimentally feasible and could be implemented using high-repetition-rate, fiber-coupled photoconductive-antenna-based THz-TDS systems. Such THz-TDS systems, commercially available products from MenloSystems, offer SNR>100dB across a bandwidth >6THz, providing an appropriate platform for experimental validation.

To address the inherent limitations of traditional THz imaging, we implemented an image reconstruction technique that leverages the complete temporal information obtained from TDS detection. This approach yields significantly better results than conventional THz imaging techniques, which typically could operate on fixed TDS delays for image reconstruction. Furthermore, the traditional method of temporal response reconstruction, which performs reconstructions at discrete delays, needs to be improved to fully characterize the optical properties of a sample, particularly when objects are composed of complex spectral responses. To illustrate this point, we integrated our approach in the spectral domain; in practice, the time-Fourier transform of the recorded TDS image is performed for each pixel, followed by solving the optimization problem defined in Eq. (3). After solving for each frequency ω, the immediate inverse time-Fourier transform of each spatial point of hyperspectral images would result in the time-resolved reconstructed images of the sample. Figure 2 demonstrates the capabilities of our full-field methodology for achieving spatiotemporal high-resolution image retrieval, particularly for phase-sensitive reconstruction of samples with complex phase profiles. As an initial demonstration, a phase screen [Fig. 2(a)] modulates the incident field, which is further recorded using TDS detection at the image plane, resulting in diffraction-limited, low-resolution spatiotemporal images [from Eq. (2)]. Figure 2(b) shows the distorted images recorded at three different times, t=−0.45, 0, and 0.63 ps, relative to our time reference of the pulse. In contrast, the image retrieval outcomes, shown in Fig. 2(c), represent the reconstructed phase images of the transmitted field at the same time values, highlighting the high fidelity of our time-resolved image reconstruction technique. After processing, the waveforms corresponding to two distinct pixels [marked by blue and red stars in Fig. 2(c)] are shown in Fig. 2(d), which are characterized by 0 and π phases of the pulses and showcase the ability to reconstruct temporal evolution of the images. Figure 2(e) illustrates the enhancement in image resolution by comparing line cuts of the phase image [dotted line in Fig. 2(a))] and the reconstructed image [dotted line in Fig. 2(c)], thereby underscoring the efficacy of our approach in improving image resolution. Interestingly, our imaging methodology can fully resolve the 150 μm slits separated by 300 μm, present in the imaging object, with marginally reduced phase contrast, which mainly depends on SNR and optical properties of the sample.

For a more challenging demonstration, our image reconstruction technique could subsequently be adapted for application to imaging in complex media. Notably, due to the larger wavelength scale of THz waves compared to optical frequencies, complex media with tunable complexity could potentially be fabricated in real experimental scenarios. Considering this, in our simulation, we decided to mimic 3D-printed, slowly varying random surfaces that operate as passive phase-modulating elements [59], keeping the close analogy to the experiments. To model the transmission properties of such random surfaces, we define a complex transfer function as follows: T˜(x,y;ω)=a0exp(iπ(nm(ω)−1)A(x,y)),(4)where nm(ω) represents frequency-dependent refractive index of the surface material, a0 is an amplitude factor, and A(x,y) is topological amplitude characterized by position-dependent relative thickness distributed over the random surface. In our random media modeling, we assumed topological amplitude as a circular random Gaussian distribution with zero mean and standard deviation σ=12. To integrate controllable features and complexity into the phase element, we applied a Gaussian filter of width Δc to the topological amplitude A(x,y). Figures 3(a) and 3(b) present imaging target random surfaces for a0=1 and nm(ω)=1.55 where we convolve white-noise distribution in A(x,y) with a Gaussian filter with a standard deviation of Δc=150μm and 600 μm, respectively, to impose a desired spatial correlation in the transfer function of random media. A notable outcome is that the overall transmission characteristics are inherently associated with the degree of complexity, specifically the spatial correlation and amplitude range of the random media model. The corresponding reconstructions of the random media at t=0ps are shown in Figs. 3(c) and 3(d), where our employed technique effectively retrieved the phase topology of the random surfaces. We also plotted a line cut to show phase variation along the x-direction [shown as dashed and solid lines in Figs. 3(a)–3(d)] in order to assess the accuracy of the reconstruction. This demonstrates the consistency in the spatial phase profile of the reconstructed image with the random phase screen. Additionally, we calculated the spatial correlation function for both the ground truth and the reconstructed random phase screens. The results are presented in Figs. 3(g) and 3(h), providing the score for the reconstruction fidelity.

We extend our approach to obtain the transmission characteristics of a hidden object placed between two random media, as illustrated in Fig. 4. We modify the forward propagation model to evaluate the wave scattering and diffraction resulting from sequential interactions with the random phase screens and the propagation of waves across free space between them. We employ the beam propagation method (BPM) [60–62] to evaluate the wave propagation through phase slices as E˜(q+1)(x,y;ω)=T˜(q+1)(x,y;ω)·℘Δz(ω)[E˜(q)(x,y;ω)],(5)where q=1,…,Q (in our case Q=3) represents the number of slices, E˜(q+1) and E˜(q) are exit field distribution in space corresponding to the (q+1)th and qth slices, respectively, T˜(q+1) is the transfer function of the (q+1)th slice, and Δz is the separation between the slices. We represent ℘Δz(ω)[−] as a mathematical operator that propagates an electric field by the distance Δz, as per the angular spectrum propagation method. We define the free-space propagator as ℘Δz(ω)[−]=Fk→x−2D{exp(−iΔz(ωc)2−kx2−ky2)·Fx→k2D{−}}, where c is the speed of light in free space, and (kx,ky) represents spatial-frequency coordinates. The BPM, functioning as a modified forward physical model for THz pulse propagation, computes the spatiotemporal wave field through sequential diffraction and multiple scattering assessments. Further additional propagation to the field exiting the final phase slice is constrained by the pupil aperture, resulting in the field on the image plane as E˜out(x,y;ω)=Op(ω)⟦℘−Δz(ω)[E˜(s)(x,y;ω]⟧,(6)where Op(ω)⟦−⟧ denotes the pupil aperture operator and ℘−Δz(ω) represents the back-propagation operator to refocus the exit wave to the physical focal plane, i.e., at the plane of the hidden imaging object, which is placed at the conjugate focus to the imaging plane of the system. We approach the image retrieval process as a least-squares optimization problem where the disparity between the recorded and estimated fields from the modified forward model can be minimized. In such cases, we modify the fitness function for three slices as $$\begin{gathered} F(T˜r)=argmin{T˜r(1),T˜r(2),T˜r(3)}∑aA∑x∑y|E˜outa(x,y;ω) \\ −ξ˜a{T˜r(1)(x,y;ω),T˜r(2)(x,y;ω),T˜r(3)(x,y;ω)}|2, \end{gathered}$$(7)where T˜r(q) represents the retrieved transfer function of the qth slice and ξ˜a{*} refers to the modified forward model operation that predicts the spatial field distribution at frequency ω on the image plane after illuminating a sample object with an incident THz pulse defined by wavevector k→a≡(kxa,kya). Further, we use an iterative approach with the accelerated distributed Nesterov gradient-descent algorithm [63] to solve the fitness function outlined in Eq. (7). As shown in Fig. 4(a), the angular projection of THz pulses interacts with an arrangement of thin slices comprising the random surfaces and hidden objects, and the complex exit field from the final slice is filtered by the pinhole, resulting in a distribution of scrambled spatiotemporal fields on the image plane. An example of a three-slice configuration is shown in Fig. 4(b), where the first and last slices are random screens containing spatial correlation of 200 μm and 800 μm, respectively, and the middle slice represents a concealed imaging object acting as a test sample. In our simulations, each phase slice is separated by Δz=12.5mm, significantly larger than the THz pulse wavelength. As a relevant case example of recovering a hidden object with complex spatial structures, initially, we placed a fast-varying phase imaging object between the two random phase screens as a test sample, as shown in Fig. 4(b). As another instance, we placed a slowly varying random phase medium containing spatially correlated (Δc=300μm) phase elements in a circular random Gaussian distribution, as shown in Fig. 4(c). The resulting reconstruction outcomes are presented in Fig. 4(d), demonstrating the successful retrieval of each slice at t=0ps, as per the time reference of the THz pulse. The reconstruction result is also shown in Fig. 4(e) for the case of slowly varying random media. Further, we analyzed the statistical characteristics of the retrieved object, as illustrated in Fig. 4(f), by calculating the probability density function for the phase elements, showing a remarkable resemblance with the ground truth object. Moreover, we evaluated the spatial correlation embedded within the topological amplitude of the random phase screen. In Fig. 4(g), we show that the retrieved phase screen exhibits a 300 μm spatial correlation, aligning with the ground truth object. It is essential to highlight that the field-sensitive image retrieval method is effective since each slice is reconstructed simultaneously with adequate reconstruction accuracy.

As a final validation, we demonstrate the capabilities of our time-resolved imaging approach for the material characterization of complex inhomogeneous objects. We show that time-resolved image reconstruction can strongly improve hyperspectral imaging, enabling us to acquire local spectral information at each spatial coordinate of the reconstructed image. Such measurements would be crucial for spectral analysis in complex semi-transparent samples, currently presented as a significantly challenging task in conventional broadband THz imaging modalities. To explore this possibility, we simulated a semi-transparent imaging target with thickness Δ=100μm composed of Teflon, Topas, and HDPE, as shown in Fig. 5(a). These materials exhibit markedly different spectral responses at THz frequencies, particularly in the 0.5–1.5 THz range. Therefore, selecting these materials allows direct differentiation of the spectral components and their spatial distribution within the sample. In terms of material composition and taking into account the material dispersion, the transfer function for the thin object can be expressed in terms of the Fabry–Perot response as T˜(x,y;ω)=t(x,y;ω)e−iϕ(x,y;ω)1−r(x,y;ω)e−2iϕ(x,y;ω),(8)where ϕ(x,y;ω)=n˜(x,y;ω)ωΔc is the function of local refractive distribution n˜(x,y;ω)=n(x,y;ω)+iρ(x,y;ω), t(x,y;ω)=4n˜(x,y;ω)(n˜(x,y;ω)+1)2 is the transmission coefficient, and r(x,y;ω)=(1−n˜(x,y;ω))2(n˜(x,y;ω)+1)2 is the reflection coefficient. The results of our simulations are shown in Fig. 5. Initially, we obtained the spatial refractive index profiles under diffraction-limited imaging conditions. As shown in Fig. 5(b), the spatial refractive index distribution of the sample at 0.5, 1, 1.5 THz retrieved by conventional imaging results in distorted images, where the differentiation between materials is obscured. Conversely, as shown in the reconstructed hyperspectral image in Fig. 5(c), our method provides high-resolution spatial features in the refractive index profile with improved material distinguishability within the sample. Such enhanced material selectivity is further validated in Fig. 5(d), where we compare the reconstructed refractive indices with the material by analyzing the spectral response within 0.5–1.5 THz and obtaining a satisfying agreement over a broad spectral range.

In this work, we have theoretically demonstrated a technique for time-resolved THz imaging based on a Fourier synthetic aperture to overcome the limitations of traditional THz imaging methods without requiring major hardware modifications. Our methodology combines a multi-angle plane wave illumination approach and a field-sensitive THz field detection enabled by state-of-the-art TDS technology. We have demonstrated that incorporating dynamic aperture synthesis within a coherent imaging framework can achieve subwavelength resolution and enable phase-sensitive reconstruction of complex, non-dispersive samples. We also have shown image retrieval even in more complex scenarios where a phase image object is hidden between two random diffusers. Our field-sensitive optimization approach is robust, effective, and compatible with experimental conditions. As a relevant case example, we demonstrated broadband imaging, successfully retrieving spatially resolved refractive index profiles for inhomogeneous and semi-transparent samples. In terms of the feasibility of experimental implementation, we propose to use a commercially available fiber-coupled photoconductive antenna THz-TDS system that would require less than 3 h for data acquisition (200×200 pixels image for 10 illumination angles) and full image reconstruction, demonstrating competitive performance while retaining accessibility and practicality compared with current state-of-the-art THz imaging techniques. Currently, the image reconstruction algorithm employs a convex optimization framework. While effective, it can be computationally intensive, particularly when processing large datasets or multiple angular illuminations. Algorithmic enhancements—such as more efficient solvers or parallel processing—could be explored to reduce computation time. The proposed method relies on multiple pulse illuminations from various angles to synthesize the aperture in Fourier space. Reducing the number of required illuminations without compromising resolution remains a challenge. Potential improvements could involve adaptive illumination strategies, which would select only the most informative angles, thereby reducing the overall number of illuminations needed for high-quality image reconstruction. Such an imaging approach offers novel insights into the optical properties and structural complexity of arbitrary samples, paving the way for new research and development opportunities across various fields, including chemical sensing, biomedical imaging, and material science using THz radiation.