# Fourier synthetic-aperture-based time-resolved terahertz imaging

Vivek Kumar, Pitambar Mukherjee, Lorenzo Valzania, Amaury Badon, Patrick Mounaix, Sylvain Gigan

► To cite this version:

Vivek Kumar, Pitambar Mukherjee, Lorenzo Valzania, Amaury Badon, Patrick Mounaix, et al. Fourier synthetic-aperture-based time-resolved terahertz imaging. Photonics research, 2025, 13 (2), pp.407-416. 10.1364/PRJ.544076. hal-04976124

HAL Id: hal-04976124

https://hal.science/hal-04976124v1

Submitted on 18 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fourier Synthetic Aperture-based Time-resolved Terahertz Imaging

Vivek Kumar<fnref n="1" />,*, Pitambar Mukherjee<fnref n="2" />, Lorenzo Valzania<fnref n="1" />, Amaury Badon<fnref n="3" />, Patrick Mounaix<fnref n="2" />, Sylvain Gigan<fnref n="1" />

**Abstract:** Terahertz microscopy has attracted attention owing to the distinctive characteristics of the THz frequency region, particularly non-ionizing photon energy, spectral fingerprint, and transparency to most nonpolar materials. However, the well-known Rayleigh diffraction limit imposed on THz waves commonly constrains the resultant imaging resolution to values beyond the millimeter scale, consequently limiting the applicability in numerous emerging applications for chemical sensing and complex media imaging. In this theoretical and numerical work, we address this challenge by introducing a new imaging approach, based on acquiring high-spatial frequencies by adapting the Fourier synthetic aperture approach to the terahertz spectral range, thus outperforming the diffraction-limited resolution. Our methodology combines multi-angle terahertz pulsed illumination with time-resolved field measurements, as enabled by the state-of-the-art time-domain spectroscopy technique. We demonstrate the potential of the approach for hyperspectral terahertz imaging of semi-transparent samples and show that the technique can reconstruct spatial and temporal features of complex inhomogeneous samples with subwavelength resolution.

**Keywords:** Terahertz Time-Domain Spectroscopy, Synthetic Aperture, Broadband Terahertz Imaging, Time-resolved Imaging.

**Abbreviations:** THz (terahertz), FOV (field-of-view), TDS (time-domain spectroscopy), SNR (signal-to-noise ratio)

## 1. Introduction:

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 systems pose 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 have 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. 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]. Similarly, ptychography imaging techniques have been introduced to accelerate high-resolution imaging

<footnote n="1">¹Laboratoire Kastler Brossel, ENS-Universite PSL, CNRS, Sorbonne Universite, College de France, 24 rue Lhomond, 75005 Paris, France.</footnote>

<footnote n="2">²IMS Laboratory, University of Bordeaux, UMR CNRS 5218, 351 Cours de la Libération Bâtiment A31, 33405 Talence, France.</footnote>

<footnote n="3">³ Laboratoire Photonique Numérique et Nanosciences (LP2N), UMR 5298, University of Bordeaux, F-33400 Talence, France.</footnote>

<footnote>*vivek.kumar@lkb.ens.fr</footnote>
without significant hardware modifications to THz imaging systems [28–30]. Quite interestingly, THz near-field imaging modalities [31–37], including tip-enhanced THz probes [35,38–40], micro-aperture THz probes [41], air-plasma dynamic aperture [42,43], 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 [44,45]. 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 detector with a finite number of pixels restricts resolution by failing to capture high-spatial frequency components. To overcome the constraints imposed by the limited space-bandwidth product in imaging systems, the synthetic aperture technique has been introduced, 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 [46,47] and subsequently adopted in optical [48–50], electron microscopy [51], sonar [52], and THz frequencies [13,53].

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 [54]. In contrast, THz time-domain spectroscopy (TDS) is an established technique, able to fully resolve electric-field oscillations in a time pulse, 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 an entirely new approach for time-resolved, coherent synthetic aperture imaging. Our methodology, which outperforms the resolution limit in conventional terahertz 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 THz beam and the sample. Consequently, our method offers hyperspectral imaging capabilities which enable the characterization of the complex, subwavelength resolution morphology of a sample across the broadband THz spectral range.

## 2. Method

### 2.1. THz-TDS imaging system based on Fourier synthetic aperture

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, $\tilde{T}(x, y; \omega)$, that we seek to
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Schematic of experimental-driven methodology. (a) Conceptual overview of methodology represents illumination of the imaging object with a sequence of tilted beams of THz pulses, and corresponding time-resolved distorted images are captured using a limited aperture THz imaging system via TDS experimental approach. The spatiotemporal fields produced by the different illumination angles are formulated and processed through a convex optimisation framework that could reconstruct high-resolution images of the sample, revealing complex details previously obscured due to the finite space-bandwidth product offered by the THz imaging system. (b) Aperture synthesis in the k-space due to different tilted beam illumination defined by wavevector $\vec{k}_a \rightarrow (k_{xa}, k_{ya})$ at frequencies 0.5, 1 and 1.5 THz (video 1). (c) Acquired distorted field amplitude associated with particular aperture shift as shown in (b). The 10 × 10 mm² object illumination area is spatially sampled at 50 µm resolution, corresponding to a number of pixels of 200 × 200.</figcaption></figure>

recover, is illuminated by a plane wave defined as $\exp(i(k_{xa}x + k_{ya}y))f(\omega)$, with the incident wave vectors $\vec{k}_a \equiv (k_{xa}, k_{ya})$ and $f(\omega)$ the THz broadband pulse spectrum, typically extending over a few THz. The transmitted field in its spatial frequency spectrum at a particular angular frequency $\omega$ can be expressed as,

$$ \tilde{E}^{+}(k_x - k_{xa}, k_y - k_{ya}; \omega) = \mathcal{F}_{x \to k}^{2D} \{ \tilde{T}(x, y; \omega) \cdot e^{i(k_{xa}x + k_{ya}y)} f(\omega) \} \quad (1) $$

where $\mathcal{F}_{x \to k}^{2D}\{-\}$ represents the 2D spatial Fourier transform and ‘.’ denotes point-wise multiplication. The transmitted field spectrum, $\tilde{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,

$$ E_{out}(x, y; t) = \mathcal{F}_{\omega \to t}^{-1D} \left[ \mathcal{F}_{k \to x}^{-2D} \{ \tilde{E}^{+}(k_x - k_{xa}, k_y - k_{ya}; \omega) \cdot O_{p}^{(\omega)}(k_x, k_y) \} \right] \quad (2) $$
where $O_p^{(\omega)}$ denotes the pupil aperture due to pinhole, $\mathcal{F}_{k \to x}^{-2D\{\cdot\}}$ represents 2D spatial inverse Fourier transform and $\mathcal{F}_{\omega \to t}^{-1D[-\cdot]}$ denotes as inverse time-Fourier transform. The overall transmitted spatiotemporal field (in Eq.(1)) at the imaging plane can be experimentally obtained using a combination of 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 $\omega$ and having a finite numerical aperture (NA), the pupil aperture is represented as a circular area formed by $k_x^2 + k_y^2 < (\frac{NA\omega}{2\pi c})^2$ in the spatial-frequency space. This means that TDS images formed at the image plane contain spatial-frequency content limited by pupil aperture. Synthetic aperture imaging exploits the fact that the $(-k_{xa}, -k_{ya})$ wavevector offset in the transmitted field images, as shown in Eq. (1,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 video 1). Fig. 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 are 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, $\tilde{E}_{out}^a(x, y; \omega)$ and the expected field $\tilde{\xi}^a$ predicted by the forward model outlined in Eq.1-2. In this scenario, we define the fitness function at $\omega$ as follows:

$$F(\tilde{T}_r) = \operatorname{argmin}_{\tilde{T}_r} \sum_{a}^{\mathcal{A}} \sum_{x} \sum_{y} |\tilde{E}_{out}^a(x, y; \omega) - \tilde{\xi}^a\{\tilde{T}_r(x, y; \omega)\}|^2 \quad (3)$$

where 'a' corresponds to the input beam having spatial wave vector $\vec{k}_a \equiv (k_{xa}, k_{ya})$, $\mathcal{A}$ is the total number of illumination angles. Here in Eq.(3), $\tilde{\xi}^a\{-\cdot\}$ 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 $\vec{k}_a$. The convex optimization problem defined in Eq. (3) is further solved by employing a Nesterov accelerated gradient-descent algorithm [55], 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. 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 100 mm² FOV is raster scanned at 50 µm resolution by the single-pixel TDS detector, mimicking experimental implementations of THz microscopy [56].

# 3. Results and discussion

## 3.1. Time-resolved field-sensitive imaging

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

<figure><img src="image_3.png" /><figcaption>**Figure 2:** Time-resolved retrieval of phase image object. (a) Phase object under investigation. (b) Measured diffraction-limited, low-resolution field images at t = -0.45, 0, and 0.63 ps obtained using conventional imaging system. (c) Reconstructed phase images at t = -0.45, 0, and 0.63 ps. (d) Temporal evolution of reconstructed THz pulses for two different pixels (red and blue dots shown in panel (c)). (e) Line-cut plot for phase imaging object (red dashed line) and reconstructed image (cyan solid line). The 10 × 10 mm² object illumination area is spatially sampled at 50 µm resolution, corresponding to a number of pixels of 200 × 200. We considered a 1 nJ THz pulse of duration 250 fs at the input with 30 dB SNR per pixel.</figcaption></figure>
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 points 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,

<figure><img src="image_4.png" /><figcaption>**Figure 3:** Time-resolved retrieval of random media. (a, b) Slowly varying random phase surfaces with spatial correlation-150 µm and 600 µm, respectively. (c, d) Retrieved spatial phase distribution of random media at t = 0 ps. (e) Line-cut in the x-direction for random phase screen (white dot line in panel a) and reconstructed phase image (red solid line in panel a). (f) line-cut in the x-direction for random phase screen (white dot line in panel b) and reconstructed phase image (green solid line in panel b). (g,h) Spatial autocorrelation function for the ground truth (black dashed) and retrieved (red and green) random phase screens. The 10 × 10 mm² object illumination area is spatially sampled at 50 µm resolution, corresponding to a number of pixels of 200 × 200.</figcaption></figure>
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 panel (c)) are shown in Fig. 2(d), which are characterized by zero and $\pi$ 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 panel (a)) and the reconstructed image (dotted line in panel (c)), thereby underscoring the efficacy of our approach in improving the 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 random media imaging. Notably, due to the larger wavelength scale of THz waves compared to optical frequencies, random media with tunable complexity could potentially be fabricated in real experimental scenarios. Considering this, in our simulation, we decided to mimic a 3D-printed, slowly varying random surfaces that operate as passive phase-modulating elements [57], keeping the close analogy to the experiments. To model the transmission properties of such a random surface, we define a complex transfer function as follows:

$$ \tilde{T}(x, y; \omega) = a_0 \exp(i\pi(n_m(\omega) - 1)A(x, y)) \quad (4) $$

where $n_m(\omega)$ represents frequency dependent refractive index of the surface material, $a_0$ 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 $\sigma = \frac{1}{2}$. To integrate controllable features and complexity into the phase element, we applied a Gaussian filter of width $\Delta c$ to the topological amplitude $A(x, y)$. Figure 3(a, b) presents imaging target random surfaces for $a_0 = 1$ and $n_m(\omega) = 1.55$ where we convolve white-noise distribution in $A(x, y)$ with a Gaussian filter with a standard deviation of $\Delta 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 within the random media model. The corresponding reconstructions of the random media at t = 0 ps are shown in Fig. 3(c, 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 Fig. 3(a-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 Fig. 3(g,h), providing the score for the reconstruction fidelity.

## 3.2. Retrieval of terahertz imaging object concealed by random media

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) [58–60] to evaluate the wave propagation through phase slices as,

$$ \tilde{E}^{(q+1)}(x, y; \omega) = \tilde{T}^{(q+1)}(x, y; \omega) \cdot \wp_{\Delta z}^{(\omega)}[\tilde{E}^{(q)}(x, y; \omega)] \quad (5) $$

where $q = 1, \dots, Q$ (in our case $Q = 3$) representing number of slices, $\tilde{E}^{(q+1)}$ and $\tilde{E}^{(q)}$ are exit field distribution in space corresponding to the $(q+1)^{th}$ and $(q)^{th}$ slices, respectively, $\tilde{T}^{q+1}$ is the transfer function of the $(q+1)^{th}$ slice and $\Delta z$ is the separation between the slices. We represent $\wp_{\Delta z}^{(\omega)}[-]$ as a mathematical operator that propagates an electric field by the distance $\Delta z$, as per the angular spectrum propagation method. We define the free space propagator as,

$$ \wp_{\Delta z}^{(\omega)}[-] = \mathcal{F}_{k \to x}^{2D} \left\{ \exp \left( -i \Delta z \sqrt{\left(\frac{\omega}{c}\right)^2 - k_x^2 - k_y^2} \right) \cdot \mathcal{F}_{x \to k}^{2D} \{-\} \right\}, \text{ where } c \text{ is the speed of light in free space, } (k_x, k_y) \text{ 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,} $$

$$ \tilde{E}_{\text{out}}(x, y; \omega) = O_p^{(\omega)} \left[ \wp_{-\Delta z}^{(\omega)}[\tilde{E}^{(s)}(x, y; \omega)] \right] \quad (6) $$

where $O_p^{(\omega)}[-]$ denotes the pupil aperture operator and $\wp_{-\Delta z}^{(\omega)}$ 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 field from the modified forward model can be minimized. In such case, we modify the fitness function for 3 slices as,

$$ \mathbb{F}(\tilde{T}_r) = \operatorname*{argmin}_{\tilde{T}_r^{(1)}, \tilde{T}_r^{(2)}, \tilde{T}_r^{(3)}} \sum_a^A \sum_x \sum_y \left| \tilde{E}_{\text{out}}^a(x, y; \omega) - \tilde{\xi}^a \left\{ \tilde{T}_r^{(1)}(x, y; \omega), \tilde{T}_r^{(2)}(x, y; \omega), \tilde{T}_r^{(3)}(x, y; \omega) \right\} \right|^2 \quad (7) $$

where $\tilde{T}_r^{(q)}$ represents the retrieved transfer function of the $(q)^{th}$ slice and $\tilde{\xi}^a\{\}$ refers to the modified forward model operation that predicts the spatial field distribution at frequency $\omega$ on the image plane after illuminating a sample object with an incident THz pulse defined by wavevector $\vec{k}_a \equiv (k_{xa}, k_{ya})$. Further, we use an iterative approach with the accelerated distributed Nesterov gradient descent algorithm [61] 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 3-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 as a test sample. In our simulations, each phase slice is separated by $\Delta z = 12.5$ mm, significantly larger than the THz pulse wavelength. The resulting reconstruction outcomes are presented in Fig. 4(c), demonstrating the successful retrieval of the hidden object at $t = 0$ ps, as per the time reference of the THz pulse. It is essential to highlight that the field-sensitive image retrieval method is effective since each slice is reconstructed simultaneously with adequate reconstruction accuracy.

Recovering a concealed object with a complex structure containing slowly varying random spatial phase patterns can be especially difficult. To explore this, we simulated our image retrieval method by placing a random phase medium containing spatially correlated ($\Delta c = 300$ µm) phase elements in circular random Gaussian distribution, as shown in Fig. 4(d),
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Retrieval of imaging object concealed by random media. (a) A schematic of an experimental-driven THz imaging system, where an imaging object is placed in the middle of the random phase screens separated by 25 mm. The wave interaction with random surfaces and hidden object leads to diffraction and scattering, which could be modelled as free-space wave propagation and wave scattering, respectively. The impinging wavefronts illuminate the same area of phase objects and provide redundant information in the resulting output waveform, further enabling the simultaneous reconstruction of the transmission properties of the random phase screen and the hidden objects. (b) Hidden object under examination with random phase surfaces containing a spatial correlation of 200 µm and 600 µm, respectively. (c) Reconstructed phase profiles of the random surfaces and the hidden imaging object at t = 0 ps. (d) Hidden random phase screen with a spatial correlation of 300 µm. (e) Retrived hidden random phase screen. (f) Histogram (Probability density function) for the phase element of ground truth (yellow shaded area) and retrieved (blue) the random phase screen. (g) Spatial autocorrelation function for the ground truth (yellow shaded area) and retrieved (blue) random phase screens. The 10 × 10 mm² object illumination area is spatially sampled at 50 µm resolution, corresponding to a number of pixels of 200 × 200.</figcaption></figure>

between the two random phase screens. The success of our approach is demonstrated in Fig. 4(e), where the retrieved phase profile of the concealed object is shown. 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.

### 3.3. Hyperspectral imaging for material characterization

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 $\Delta = 100 \ \mu 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 Fabry-Perot response as,

$$ \tilde{T}(x, y; \omega) = \frac{t(x, y; \omega)e^{-i\phi(x, y; \omega)}}{1 - r(x, y; \omega)e^{-2i\phi(x, y; \omega)}} \quad (8) $$

where $\phi(x, y; \omega) = \frac{\tilde{n}(x, y; \omega)\omega\Delta}{c}$ is the function of local refractive distribution $\tilde{n}(x, y; \omega) = n(x, y; \omega) + i\rho(x, y; \omega)$, $t(x, y; \omega) = \frac{4\tilde{n}(x, y; \omega)}{[\tilde{n}(x, y; \omega)+1]^2}$ is transmission coefficient and $r(x, y; \omega) =$

<figure><img src="image_6.png" /><figcaption>Figure 5: Hyperspectral imaging for material characterisation. (a) The schematic of the 100 µm thick imaging object composed of Teflon, Topas and HDPE, showing the spatial distribution of the refractive index at 1 THz. (b) Distorted spatial variation in the refractive index profile at 0.5, 1 and 1.5 THz recorded using a diffraction-limited imaging system. (c) Reconstructed spatial refractive index profiles at 0.5, 1 and 1.5 THz illustrating the high-resolution material selectivity. (d) Variation in the refractive index of Teflon, Topas, and HDPE as a function of frequency representing comparison in the reconstructed refractive indices (solid lines) and material refractive indices (dashed lines). The 10 × 10 mm² object illumination area is spatially sampled at 50 µm resolution, corresponding to a number of pixels of 200 × 200.</figcaption></figure>
$[1-\tilde{n}(x,y;\omega)]^2 / [\tilde{n}(x,y;\omega)+1]^2$ is 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 spectral response within 0.5-1.5 THz and obtaining a satisfying agreement over a broad spectral range.

## 4. Conclusion

In this work, we have theoretically demonstrated a technique for time-resolved THz imaging based on 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 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. 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 debate. 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, opening new avenues for research and development in various fields, including chemical sensing, biomedical imaging, and material science using THz radiation.

**Funding:** This project received funding from HYPSTER project under the ANR within the France 2030 framework: ANR-22-CE42-0005 and ANR 22-PEEL-0003 Comptera.

## Supplementary materials

The Supporting Information is available free of charge at [LINK].

The datasets for all figures are freely accessible at: [LINK-figshare]

Additional videos for aperture synthesis in the Fourier space and corresponding image formation (video 1).

## Author Contributions

All authors were engaged in the general discussion regarding the basic science of the paper. V.K. conceived the idea, performed the calculations reported and drafted the paper. All authors contributed to the general understanding of the results. P.Mo. and S.G. supervised the general research activities.
**Disclosures:** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are publicly available.
Additional data can be requested from the corresponding author.

**References**

1. B. Recur, H. Balacey, J. B. Sleiman, J. B. Perraud, J.-P. Guillet, A. Kingston, and P. Mounaix, "Ordered subsets convex algorithm for 3D terahertz transmission tomography," Opt. Express, OE 22, 23299-23309 (2014).

2. J. B. Sleiman, B. Bousquet, N. Palka, and P. Mounaix, "Quantitative Analysis of Hexahydro-1,3,5-trinitro-1,3,5, Triazine/Pentaerythritol Tetranitrate (RDX-PETN) Mixtures by Terahertz Time Domain Spectroscopy," Appl Spectrosc 69, 1464-1471 (2015).

3. P. Hillger, R. Jain, J. Grzyb, W. Förster, B. Heinemann, G. MacGrogan, P. Mounaix, T. Zimmer, and U. R. Pfeiffer, "A 128-Pixel System-on-a-Chip for Real-Time Super-Resolution Terahertz Near-Field Imaging," IEEE Journal of Solid-State Circuits 53, 3599-3612 (2018).

4. K. Okada, K. Serita, Q. Cassar, H. Murakami, G. MacGrogan, J.-P. Guillet, P. Mounaix, and M. Tonouchi, "Terahertz near-field microscopy of ductal carcinoma in situ (DCIS) of the breast," J. Phys. Photonics 2, 044008 (2020).

5. F. Fauquet, F. Galluzzi, P. F. Taday, R. Chapoulie, A. Mounier, A. Ben Amara, and P. Mounaix, "Terahertz time-domain spectro-imaging and hyperspectral imagery to investigate a historical Longwy glazed ceramic," Sci Rep 14, 19248 (2024).

6. A. W. Lohmann, R. G. Dorsch, D. Mendlovic, Z. Zalevsky, and C. Ferreira, "Space-bandwidth product of optical signals and systems," J. Opt. Soc. Am. A, JOSAA 13, 470-473 (1996).

7. S.-H. Ding, Q. Li, Y.-D. Li, and Q. Wang, "Continuous-wave terahertz digital holography by use of a pyroelectric array camera," Opt. Lett., OL 36, 1993-1995 (2011).

8. P. Zolliker and E. Hack, "THz holography in reflection using a high resolution microbolometer array," Opt. Express, OE 23, 10957-10967 (2015).

9. M. Locatelli, M. Ravaro, S. Bartalini, L. Consolino, M. S. Vitiello, R. Cicchi, F. Pavone, and P. De Natale, "Real-time terahertz digital holography with a quantum cascade laser," Sci Rep 5, 13566 (2015).

10. M. Yamagiwa, T. Ogawa, T. Minamikawa, D. G. Abdelsalam, K. Okabe, N. Tsurumachi, Y. Mizutani, T. Iwata, H. Yamamoto, and T. Yasui, "Real-Time Amplitude and Phase Imaging of Optically Opaque Objects by Combining Full-Field Off-Axis Terahertz Digital Holography with Angular Spectrum Reconstruction," J Infrared Milli Terahz Waves 39, 561-572 (2018).

11. M. Humphreys, J. P. Grant, I. Escorcia-Carranza, C. Accarino, M. Kenney, Y. D. Shah, K. G. Rew, and D. R. S. Cumming, "Video-rate terahertz digital holographic imaging system," Opt. Express, OE 26, 25805-25813 (2018).

12. K. Xue, Q. Li, Y.-D. Li, and Q. Wang, "Continuous-wave terahertz in-line digital holography," Opt. Lett., OL 37, 3228-3230 (2012).

13. H. Huang, L. Rong, D. Wang, W. Li, Q. Deng, B. Li, Y. Wang, Z. Zhan, X. Wang, and W. Wu, "Synthetic aperture in terahertz in-line digital holography for resolution enhancement," Appl. Opt., AO 55, A43-A48 (2016).

14. Z. Li, L. Li, Y. Qin, G. Li, D. Wang, and X. Zhou, "Resolution and quality enhancement in terahertz in-line holography by sub-pixel sampling with double-distance reconstruction," Opt. Express, OE 24, 21134-21146 (2016).

15. H. Huang, D. Wang, L. Rong, S. Panezai, D. Zhang, P. Qiu, L. Gao, H. Gao, H. Zheng, and Z. Zheng, "Continuous-wave off-axis and in-line terahertz digital holography with
phase unwrapping and phase autofocusing," *Optics Communications* **426**, 612-622 (2018).

16. Z. Li, Q. Yan, Y. Qin, W. Kong, G. Li, M. Zou, D. Wang, Z. You, and X. Zhou, "Sparsity-based continuous wave terahertz lens-free on-chip holography with sub-wavelength resolution," *Opt. Express*, OE **27**, 702-713 (2019).

17. L. Rong, T. Latychevskia, C. Chen, D. Wang, Z. Yu, X. Zhou, Z. Li, H. Huang, Y. Wang, and Z. Zhou, "Terahertz in-line digital holography of human hepatocellular carcinoma tissue," *Sci Rep* **5**, 8445 (2015).

18. L. Olivieri, J. S. Totero Gongora, A. Pasquazi, and M. Peccianti, "Time-Resolved Nonlinear Ghost Imaging," *ACS Photonics* **5**, 3379-3388 (2018).

19. J. Zhao, Y. E. K. Williams, X.-C. Zhang, and R. W. Boyd, "Spatial sampling of terahertz fields with sub-wavelength accuracy via probe-beam encoding," *Light: Science & Applications* **8**, 55 (2019).

20. J. S. Totero Gongora, L. Olivieri, L. Peters, J. Tunesi, V. Cecconi, A. Cutrona, R. Tucker, V. Kumar, A. Pasquazi, and M. Peccianti, "Route to Intelligent Imaging Reconstruction via Terahertz Nonlinear Ghost Imaging," *Micromachines* **11**, 521 (2020).

21. M. P. L. Olivieri, J. S. T. Gongora, L. Peters, V. Cecconi, A. Cutrona, J. Tunesi, R. Tucker, A. Pasquazi, "Hyperspectral terahertz microscopy via nonlinear ghost imaging," *Optica* **7**, 186 (2020).

22. R. I. Stantchev, X. Yu, T. Blu, and E. Pickwell-MacPherson, "Real-time terahertz imaging with a single-pixel detector," *Nat Commun* **11**, 2535 (2020).

23. V. Kumar, V. Cecconi, L. Peters, J. Bertolotti, A. Pasquazi, J. S. Totero Gongora, and M. Peccianti, "Deterministic Terahertz Wave Control in Scattering Media," *ACS Photonics* **9**, 2634-2642 (2022).

24. V. Cecconi, V. Kumar, A. Pasquazi, J. S. Totero Gongora, and M. Peccianti, "Nonlinear field-control of terahertz waves in random media for spatiotemporal focusing," *Open Res Europe* **2**, 32 (2022).

25. L. Olivieri, L. Peters, V. Cecconi, A. Cutrona, M. Rowley, J. S. Totero Gongora, A. Pasquazi, and M. Peccianti, "Terahertz Nonlinear Ghost Imaging via Plane Decomposition: Toward Near-Field Micro-Volumetry," *ACS Photonics* **10**, 1726-1734 (2023).

26. V. Cecconi, V. Kumar, J. Bertolotti, L. Peters, A. Cutrona, L. Olivieri, A. Pasquazi, J. S. Totero Gongora, and M. Peccianti, "Terahertz Spatiotemporal Wave Synthesis in Random Systems," *ACS Photonics* (2024).

27. X. Li, J. Li, Y. Li, A. Ozcan, and M. Jarrahi, "High-throughput terahertz imaging: progress and challenges," *Light Sci Appl* **12**, 233 (2023).

28. L. Valzania, T. Feurer, P. Zolliker, and E. Hack, "Terahertz ptychography," *Opt. Lett.*, OL **43**, 543-546 (2018).

29. L. Valzania, P. Zolliker, and E. Hack, "Coherent reconstruction of a textile and a hidden object with terahertz radiation," *Optica, OPTICA* **6**, 518-523 (2019).

30. L. Rong, C. Tang, D. Wang, B. Li, F. Tan, Y. Wang, and X. Shi, "Probe position correction based on overlapped object wavefront cross-correlation for continuous-wave terahertz ptychography," *Opt. Express*, OE **27**, 938-950 (2019).

31. S. Hunsche, "THz near-field imaging," *Optics Communications* **150**, 22-26 (1998).

32. Q. Chen, Z. Jiang, G. X. Xu, and X. C. Zhang, "Near-field terahertz imaging with a dynamic aperture.," *Optics letters* **25**, 1122-1124 (2000).

33. N. C. J. van der Valk and P. C. M. Planken, "Electro-optic detection of subwavelength terahertz spot sizes in the near field of a metal tip," *Appl. Phys. Lett.* **81**, 1558-1560 (2002).

34. H.-T. Chen, R. Kersting, and G. C. Cho, "Terahertz imaging with nanometer resolution," *Applied Physics Letters* **83**, 3009-3009 (2003).
35. A. J. Huber, F. Keilmann, J. Wittborn, J. Aizpurua, and R. Hillenbrand, "Terahertz near-field nanoscopy of mobile carriers in single semiconductor nanodevices," Nano Letters **8**, 3766-3770 (2008).

36. X. Chen, X. Liu, X. Guo, S. Chen, H. Hu, E. Nikulina, X. Ye, Z. Yao, H. A. Bechtel, M. C. Martin, G. L. Carr, Q. Dai, S. Zhuang, Q. Hu, Y. Zhu, R. Hillenbrand, M. Liu, and G. You, "THz Near-Field Imaging of Extreme Subwavelength Metal Structures," ACS Photonics **7**, 687-694 (2020).

37. L. L. Hale, T. Siday, and O. Mitrofanov, "Near-field imaging and spectroscopy of terahertz resonators and metasurfaces [Invited]," Opt. Mater. Express, OME **13**, 3068-3086 (2023).

38. T. L. Cocker, V. Jelic, M. Gupta, S. J. Molesky, J. A. J. Burgess, G. D. L. Reyes, L. V. Titova, Y. Y. Tsui, M. R. Freeman, and F. A. Hegmann, "An ultrafast terahertz scanning tunneling microscope," Nature Photon **7**, 620-625 (2013).

39. O. Mitrofanov, L. Viti, E. Dardanis, M. C. Giordano, D. Ercolani, A. Politano, L. Sorba, and M. S. Vitiello, "Near-field terahertz probes with room-temperature nanodetectors for subwavelength resolution imaging," Sci Rep **7**, 44240 (2017).

40. X. Guo, K. Bertling, B. C. Donose, M. Brünig, A. Cernescu, A. A. Govyadinov, and A. D. Rakić, "Terahertz nanoscopy: Advances, challenges, and the road ahead," Applied Physics Reviews **11**, 021306 (2024).

41. Y. Kawano and K. Ishibashi, "An on-chip near-field terahertz probe and detector," Nature Photonics **2**, 618-621 (2008).

42. X. Wang, J. Ye, W. Sun, P. Han, L. Hou, and Y. Zhang, "Terahertz near-field microscopy based on an air-plasma dynamic aperture," Light Sci Appl **11**, 129 (2022).

43. D. Headland and W. Withayachumnankul, "In the shadow of the laser phantom needle cross: dynamic air-plasma aperture sheds light on terahertz microscopy," Light Sci Appl **11**, 147 (2022).

44. K. Khare, *Fourier Optics and Computational Imaging* (Wiley, John Wiley & Sons Ltd, 2016).

45. J. W. Goodman, *Statistical Optics* (John Wiley & Sons, 2015).

46. M. Ryle and A. Hewish, "The Synthesis of Large Radio Telescopes," Monthly Notices of the Royal Astronomical Society **120**, 220-230 (1960).

47. H. Touma, "Synthetic Aperture Technic in Astronomy Using Slit Aperture Telescope," Astrophysics and Space Science **258**, 47-53 (1997).

48. F. L. Clerc, M. Gross, and L. Collot, "Synthetic-aperture experiment in the visible with on-axis digital heterodyne holography," Opt. Lett., OL **26**, 1550-1552 (2001).

49. T. Gutzler, T. R. Hillman, S. A. Alexandrov, and D. D. Sampson, "Coherent aperture-synthesis, wide-field, high-resolution holographic microscopy of biological tissue," Opt. Lett., OL **35**, 1136-1138 (2010).

50. W. Luo, A. Greenbaum, Y. Zhang, and A. Ozcan, "Synthetic aperture-based on-chip microscopy," Light Sci Appl **4**, e261-e261 (2015).

51. T. Hidaka, "Method for synthetic-aperture electron microscopy," J. Opt. Soc. Am., JOSA **66**, 147-150 (1976).

52. R. E. Hansen, "Introduction to Synthetic Aperture Sonar," in *Sonar Systems* (IntechOpen, 2011).

53. H. S. Lui, T. Taimre, K. Bertling, Y. L. Lim, P. Dean, S. P. Khanna, M. Lachab, A. Valavanis, D. Indjin, E. H. Linfield, A. G. Davies, and A. D. Rakić, "Terahertz inverse synthetic aperture radar imaging using self-mixing interferometry with a quantum cascade laser," Opt. Lett., OL **39**, 2629-2632 (2014).

54. N. V. Petrov, J.-B. Perraud, A. Chopard, J.-P. Guillet, O. A. Smolyanskaya, and P. Mounaix, "Terahertz phase retrieval imaging in reflection," Opt. Lett., OL **45**, 4168-4171 (2020).
55. G. Qu and N. Li, "Accelerated Distributed Nesterov Gradient Descent," IEEE Transactions on Automatic Control **65**, 2566-2581 (2020).

56. M. Naftaly, "Metrology Issues and Solutions in THz Time-Domain Spectroscopy: Noise, Errors, Calibration," IEEE Sensors Journal **13**, 8-17 (2013).

57. C. H. Brodie, I. Spotts, H. Reguigui, C. A. Leclerc, M. E. Mitchell, J. F. Holzman, and C. M. Collier, "Comprehensive study of 3D printing materials over the terahertz regime: absorption coefficient and refractive index characterizations," Opt. Mater. Express, OME **12**, 3379-3402 (2022).

58. A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, "Ptychographic transmission microscopy in three dimensions using a multi-slice approach," J. Opt. Soc. Am. A, JOSAA **29**, 1606-1614 (2012).

59. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, "Optical Tomographic Image Reconstruction Based on Beam Propagation and Sparse Regularization," IEEE Transactions on Computational Imaging **2**, 59-70 (2016).

60. X. Ma, W. Xiao, and F. Pan, "Optical tomographic reconstruction based on multi-slice wave propagation method," Opt. Express, OE **25**, 22595-22607 (2017).

61. G. Qu and N. Li, "Accelerated distributed Nesterov Gradient Descent for convex and smooth functions," in *2017 IEEE 56th Annual Conference on Decision and Control (CDC)* (2017), pp. 2260-2267.