Performance Evaluation of a THz Pulsed Imaging System: Point Spread Function, Broadband THz Beam Visualization and Image Reconstruction

: Terahertz (THz) technology is a promising research ﬁeld for various applications in basic science and technology. In particular, THz imaging is a new ﬁeld in imaging science, where theories, mathematical models and techniques for describing and assessing THz images have not completely matured yet. In this work, we investigate the performances of a broadband pulsed THz imaging system (0.2–2.5 THz). We characterize our broadband THz beam, emitted from a photoconductive antenna (PCA), and estimate its point spread function (PSF) and the corresponding spatial resolution. We provide the ﬁrst, to our knowledge, 3D beam proﬁle of THz radiation emitted from a PCA, along its propagation axis, without the using of THz cameras or proﬁlers, showing the beam spatial intensity distribution. Finally, we evaluate the THz image formation on a test-sample composed by a regular linen natural pattern.

In particular, THz time-domain spectroscopy (THz-TDS) is based on the synchronous and coherent detection of a THz pulse: both its amplitude and phase are recorded as a function of time, allowing the measurement of sample optical properties, such as the refractive index and the absorption coefficient, without the need of Kramers-Kronig relations.
Actually, THz radiation can be produced and detected from a variety of designed THz emitters and receivers, both in continuous wave (CW) and in pulsed regime. Among

Experimental Setup
Our home-made TPI system is based on a pump-probe configuration in transmission mode that produces broadband THz radiation in the range 0.2-3 THz. The experimental layout is illustrated in Figure 1.
Two twin PCAs (G10620-11, Hamamatsu) are used for both THz generation and detection. The system is powered by a mode-locked femtosecond laser (FemtoFiberNIRpro, Toptica) at 780 nm, with 100 fs pulse duration and 80 MHz repetition rate. A 50:50 beamsplitter splits the laser beam into two parts: the pump and the probe beams (each 15 mW in power). Dielectric mirrors (M, BB1-E03, Thorlabs) convey the two laser beams towards emitter and receiver PCAs. The pump beam is modulated by the THz emitter bias voltage (−7/7 V, 1 kHz) and then focused on the THz PCA emitter. The generated THz radiation is then collimated and focused onto a target-sample by using two polymethylpentene plano-convex lenses (TPX, MenloSystems, f = 50 mm). After the sample, the transmitted THz beam is recollimated and refocused onto the THz PCA detector using two TPX lenses (f = 50 mm). Simultaneously, the probe beam is used to gate the THz PCA detector. A motorized delay line (DL, DDSM100/M Thorlabs) is used to temporally and coherently sample the THz electric field, and thus both the amplitude and phase are measured (See Section S1 in Supplementary Materials). The signal from PCA detector is filtered and detected by a lock-in amplifier (LIA, SR530, Stanford). The analog output of LIA is collected and digitalized by a dedicated National Instrument acquisition card at 16-bit (DAQ, NI, 6361-BNC connector) and then transferred to a computer for data collection, visualization and off-line analysis. In order to collect THz images, the sample is aligned perpendicularly to the propagating THz radiation and can be moved in the THz beam, by a 3D imaging module. It is formed of a motorized 3D linear translation stage (PT1/M, Thorlabs with 25 mm travel range and 10 µm accuracy and servo actuators Z825B, Thorlabs with <5 µm repeatability) that moves the specimen along the THz beam propagation axis (z-axis, 1D motion) and perpendicularly to it (x,y-axes, 2D motion). The object surface is sampled on a discrete grid composed of square pixels and an acquisition for each pixel is made. Therefore, the measurement is performed with a step-by-step sampling and it can be acquired in fixed DL or spectral modes. For fixed DL images, the DL is set to a fixed position, for example at the location that corresponds to the maximum of the THz electric field E MAX . Thus, only E MAX is acquired in each pixel. For spectral images, instead, the information of a complete DL temporal scan is recorded, in each pixel. It allows the reconstruction of the whole THz pulse. By means of Fast Fourier Transform (FFT), images at given frequencies can be extracted from spectral images. Two twin PCAs (G10620-11, Hamamatsu) are used for both THz generation and detection. The system is powered by a mode-locked femtosecond laser (FemtoFiberNIRpro, Toptica) at 780 nm, with 100 fs pulse duration and 80 MHz repetition rate. A 50:50 beamsplitter splits the laser beam into two parts: the pump and the probe beams (each 15 mW in power). Dielectric mirrors (M, BB1-E03, Thorlabs) convey the two laser beams towards emitter and receiver PCAs. The pump beam is modulated by the THz emitter bias voltage (-7/7 V, 1 kHz) and then focused on the THz PCA emitter. The generated THz radiation is then collimated and focused onto a target-sample by using two polymethylpentene planoconvex lenses (TPX, MenloSystems, f = 50 mm). After the sample, the transmitted THz beam is recollimated and refocused onto the THz PCA detector using two TPX lenses (f = 50 mm). Simultaneously, the probe beam is used to gate the THz PCA detector. A motorized delay line (DL, DDSM100/M Thorlabs) is used to temporally and coherently sample the THz electric field, and thus both the amplitude and phase are measured (See Section S1 in Supplementary materials). The signal from PCA detector is filtered and detected by a lock-in amplifier (LIA, SR530, Stanford). The analog output of LIA is collected and digitalized by a dedicated National Instrument acquisition card at 16-bit (DAQ, NI, 6361-BNC connector) and then transferred to a computer for data collection, visualization and

Theory
In this section, we briefly recall the basic concepts of spatial resolution and point spread function (PSF), and later, we report the main quantities for a Gaussian beam. These calculations will be the fundamentals to understand the experimental results found in our work.

Spatial Resolution
The theoretical spatial resolution R T at a certain frequency ( f ) or wavelength (λ) is generally limited by the numerical aperture (N A) of the imaging system and can be defined by the Rayleigh limit [41]  where λ = c/ f is the wavelength and c the speed of light in vacuum. As a consequence, R T ( f ) is strongly frequency-dependent and small objects or any spatial features are easier resolvable at high frequencies. Therefore, in a broadband THz imaging system it is crucial to calculate the spatial resolution for each frequency of interest.

Point Spread Function
PSF is defined as the imaging system response to an ideal, point-like source and constitutes the basic unit forming the image [27,28,42,43] Thus, it represents the optical system's fidelity in reproducing objects through their images. In a perfect imaging system, the intensity coming from a point source (in the object plane) would be concentrated at a single point (in the image plane), the ideal image point. In real systems, however, optical imperfections result in an intensity distribution broadening around the ideal image point, producing blurred images. The PSF provides an estimation of this blurring effect [42,43]. Considering a linear system, the output formed by a sum of inputs is equal to the sum of the outputs corresponding to each input acting separately. Mathematically, if an object o(x, y) is imaged with an imaging system having a response PSF(x, y), the observed image i(x, y) is given by: where x and y are plane coordinates of the imaging system, perpendicular to the propagation beam axis and * denotes convolution operation [27]. As i(x, y) is the acquired image, Equation (2) can be inversely used in order to calculate o(x, y), knowing the system PSF(x, y). PSF of the imaging system is therefore used for improving image quality. Moreover, in a typical TPI system, the THz pulse is broadband, and thus it cannot be treated as a monochromatic beam. Therefore, the PSF is a frequency-dependent quantity: for including the full spectrum, the PSF is reconstructed by superposition of monochromatic beams in the frequency band. PSF can be evaluated by imaging objects such as tiny reflecting balls [28] and subwavelength holes for reflection and transmission setup geometries, respectively.
According to Figure 2, the Gaussian beam parameters can be defined and linked to each other. W 0 is the beam waist size at the focal plane (z = 0). The beam waist W(z), along the z direction, enlarges and can be expressed by the following equation: where Z R is called Rayleigh distance and represents the distance along the z axis where W(Z R ) = √ 2W 0 . It is defined as follows where f is the frequency and c the speed of light in vacuum. Note that, as shown in Figure 2, the waist sizes, both in the focal region (W 0 ) and outside (W(z)), refer to the Gaussian width r associated to an intensity reduction of a factor 1/e 2 . Moreover, the beam divergence angle θ can be defined as the ratio between W 0 and Z R : where is the frequency and the speed of light in vacuum. Note that, as shown in Figure 2, the waist sizes, both in the focal region ( ) and outside ( ( )), refer to the Gaussian width associated to an intensity reduction of a factor 1 ⁄ . Moreover, the beam divergence angle can be defined as the ratio between and :

PSF and Spatial Resolution
As briefly recalled in Section 3.2, the PSF of an imaging system can be evaluated from sub-wavelength pinhole images [27,28,42,48]. We used a 200 ± 6 µm diameter metallic pinhole to simulate a point-like source. We put it in the focal region of our beam ( = 0, with respect to Figure 2). Note that TPX lenses used for focusing THz radiation have constant refractive index in our spectral range (n = 1.456 ± 0.001) [49]. Therefore, chromatic aberrations are negligible in our setup. Figure 3a shows four spectral images at 0.6, 1.0, 1.6 and 2.5 THz.

PSF and Spatial Resolution
As briefly recalled in Section 3.2, the PSF of an imaging system can be evaluated from sub-wavelength pinhole images [27,28,42,48]. We used a 200 ± 6 µm diameter metallic pinhole to simulate a point-like source. We put it in the focal region of our beam (z = 0, with respect to Figure 2). Note that TPX lenses used for focusing THz radiation have constant refractive index in our spectral range (n = 1.456 ± 0.001) [49]. Therefore, chromatic aberrations are negligible in our setup. Figure 3a shows four spectral images at 0.6, 1.0, 1.6 and 2.5 THz. Point-like source is seen very differently depending on what frequency is selected. In order to evaluate the PSF, the spot intensity profile is plotted as a function of distance (black points in Figure 3b, data shown for 1.0 THz image). The Full Width at Half Maximum (FWHM) has been calculated from the Gaussian fit for each image in Figure 3a. With this method, we estimated how the point-like source is seen from the imaging system. This coincides with the theoretical spatial resolution at a given frequency estimated for our imaging system. Figure 3b reports the experimental data (black) and the Gaussian fit (green) performed to extract the FWHM, at 1.0 THz. Gray level refers to the digital information contained in each pixel. It is a discrete value between 0 and 65565, because images Point-like source is seen very differently depending on what frequency is selected. In order to evaluate the PSF, the spot intensity profile is plotted as a function of distance (black points in Figure 3b, data shown for 1.0 THz image). The Full Width at Half Maximum (FWHM) has been calculated from the Gaussian fit for each image in Figure 3a. With this method, we estimated how the point-like source is seen from the imaging system. This coincides with the theoretical spatial resolution at a given frequency estimated for our imaging system. Figure 3b reports the experimental data (black) and the Gaussian fit (green) performed to extract the FWHM, at 1.0 THz. Gray level refers to the digital information contained in each pixel. It is a discrete value between 0 and 65565, because images have been acquired with a 16-bit data acquisition card (see Section 2 and Section S2 in Supplementary Materials).
Note that the Gaussian fit has been performed using the Gaussian function f (x) where µ and σ are the Gaussian expected value and variance, respectively used as fit parameters. According to [47], the FWHM is: The Table 1 below lists the FWHMs extracted from Figure 3a at various frequencies. As can be seen in Table 1, the FWHM decreases as the frequency increases. We adopted the procedure explained above for all the images at fixed frequencies in the range 0.5-2.5 THz. In Figure 4, we report the extracted FWHM, i.e., the spatial resolution, at different frequencies (black points) and the green curve fitting (Equation (1)).  As can be seen in Table 1, the FWHM decreases as the frequency increases. We adopted the procedure explained above for all the images at fixed frequencies in the range 0.5-2.5 THz. In Figure 4, we report the extracted FWHM, i.e., the spatial resolution, at different frequencies (black points) and the green curve fitting (Equation (1)).  (6) and (7), at different frequencies. Black: experimental data. Green: our imaging system spatial resolution curve. It was obtained by fitting the data with the Rayleigh limit Equation (1).
As expected, the measured resolution data points show a hyperbolic trend as the frequency increases and it is in total agreement with Equation (1). The spatial resolution is higher as the frequency increases. Note that we selected the frequency range 0.5-2.5 THz because for low frequencies (f < 0.5 THz) the diffractive effects caused by the pinhole were too strong and no luminous spot was clearly identifiable; for high frequencies (f > 2.5 THz) the signal-to-noise ratio (SNR) [50,51] was too low and no detectable intensity was clearly  (6) and (7), at different frequencies. Black: experimental data. Green: our imaging system spatial resolution curve. It was obtained by fitting the data with the Rayleigh limit Equation (1).
As expected, the measured resolution data points show a hyperbolic trend as the frequency increases and it is in total agreement with Equation (1). The spatial resolution is higher as the frequency increases. Note that we selected the frequency range 0.5-2.5 THz because for low frequencies (f < 0.5 THz) the diffractive effects caused by the pinhole were too strong and no luminous spot was clearly identifiable; for high frequencies (f > 2.5 THz) the signal-to-noise ratio (SNR) [50,51] was too low and no detectable intensity was clearly measurable.
By fitting the data with Equation (1), the NA experimental value for our broadband imaging system has been extracted and estimated (NA = 0.27 ± 0.01).

THz 3D Beam Image
In order to realize a 3D visualization of THz beam emitted from the PCA, we spatially selected the beam and acquired fixed DL images at different z positions. We used a 1.20 ± 0.05 mm diameter metallic pinhole, moved along x, y and z directions. Each image (60 px × 60 px, 15 mm × 15 mm) was acquired in-plane (x,y) at fixed z positions, collecting a z stack-image. The z step was set to 1 mm. The spot size for each image was extracted by a Gaussian fit. According to the following equation, the spot radius r s is defined as In Figure 5, we report the procedure for the image located at z = −3 mm (3 mm before the THz beam focal plane), as an example.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 14 In Figure 5, we report the procedure for the image located at = 3 mm (3 mm before the THz beam focal plane), as an example.  Figure 5b shows its intensity profile as a function of distance (black dots) and the Gaussian fit (green curve). The 3D Volume Viewer tool by ImageJ was used to convert the z-stack image into a 3D volume, in Figure 6.   Figure 5b shows its intensity profile as a function of distance (black dots) and the Gaussian fit (green curve). The 3D Volume Viewer tool by ImageJ was used to convert the z-stack image into a 3D volume, in Figure 6.
As expected, the maximum intensity locates where the THz beam has its smallest diameter, i.e., the focal region. Figure 6 is formed of a stack of fixed DL images, therefore no information about THz frequencies spatial distribution is known. The 3D image in Figure 6 is a result of all frequency contributions, each one with its own weight (See Figure S2 in Supplementary Materials for THz power spectrum vs. frequency).
The knife edge method was adopted for reconstructing the 1D THz beam profile [46,52,53] and for comparing the results with the pinhole method. A sharp metallic edge was placed vertically within the THz beam and was moved along x and z directions (the edge position along y was kept fixed for all the knife edge measurements). Each fixed DL image was acquired by varying x position, at different z-fixed positions. The z step was set to 1 mm. mm (3 mm before the focal point). (b) Intensity plot as a function of distance. The scale bar is 1 mm. Figure 5a represents a crop of the acquired image at = 3 mm in grey look-up table (LUT). Figure 5b shows its intensity profile as a function of distance (black dots) and the Gaussian fit (green curve). The 3D Volume Viewer tool by ImageJ was used to convert the z-stack image into a 3D volume, in Figure 6. As expected, the maximum intensity locates where the THz beam has its smallest diameter, i.e., the focal region. Figure 6 is formed of a stack of fixed DL images, therefore no information about THz frequencies spatial distribution is known. The 3D image in Figure 6 is a result of all frequency contributions, each one with its own weight (See Figure  S2 in Supplementary materials for THz power spectrum vs. frequency).
The knife edge method was adopted for reconstructing the 1D THz beam profile [46,52,53] and for comparing the results with the pinhole method. A sharp metallic edge In general, a sharp edge moving in a light beam having electric field E, can be modeled by a function containing the error function er f (x) [52]: where E(x, z) is the THz electric field as a function of distances along x and z directions, E MAX is the maximum THz electric field, x 50% the position along x corresponding to E MAX /2 and W(z) is the THz beam waist as a function of distance z, along the propagation axis. This equation weights all the frequency contributions because no frequency is selected in fixed DL images. By fitting the data points with Equation (9), the green curve is obtained and displayed in Figure 7.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 14 was placed vertically within the THz beam and was moved along and directions (the edge position along was kept fixed for all the knife edge measurements). Each fixed DL image was acquired by varying position, at different -fixed positions. The step was set to 1 mm. In general, a sharp edge moving in a light beam having electric field , can be modeled by a function containing the error function ( ) [52]: where ( , ) is the THz electric field as a function of distances along and directions, is the maximum THz electric field, % the position along corresponding to 2 ⁄ and ( ) is the THz beam waist as a function of distance , along the propagation axis. This equation weights all the frequency contributions because no frequency is selected in fixed DL images. By fitting the data points with Equation (9), the green curve is obtained and displayed in Figure 7. This was performed for each image at each position. The beam waist size ( ), at different positions, have been extracted from the fitting parameters, using Equation (9). Figure 8 shows the results for the two methods adopted for the beam profile reconstruction. This was performed for each image at each z position. The beam waist size W(z), at different z positions, have been extracted from the fitting parameters, using Equation (9). Figure 8 shows the results for the two methods adopted for the beam profile reconstruction.
points as a function of distance along the direction (perpendicular to the beam propagation axis). Green: fitting curve using Equation (9). This was performed for each image at each position. The beam waist size ( ), at different positions, have been extracted from the fitting parameters, using Equation (9). Figure 8 shows the results for the two methods adopted for the beam profile reconstruction. Figure 8. Comparison between the two methods adopted for THz beam profiling (hole and knife edge methods) and curve fit. In blue: hole, in orange: knife, in green: fitting curve using Equation (3). In particular, the quantity 2W(z) is plotted in Figure 8, and therefore the THz beam 1D projection is shown 1.2 cm before and after the focal region. The blue circles represent the pinhole method, the oranges show the knife edge method: they are well in agreement. By fitting the data points with Equation (3), the relevant parameters for a Gaussian beam, defined in Section 3.3, have been found. They are listed in Table 2.

Sample Test: Linen
After studying the spatial distribution properties and the capability to resolve objects of our TPI system, we used this experimental apparatus on a sample test: a natural linen pattern. This is a quite regular pattern with horizontal and vertical linen fibers, separated by air holes (see Figure 9f). Each vertical fiber width is~0.40 ± 0.05 mm (~20 fibers/cm) and each horizontal fiber~0.15 ± 0.05 mm (~17 fibers/cm). Holes sizes vary from~50 µm to~600 µm.
Linen spectral images have been acquired. We collected transmittance images at different frequencies. Each pixel contains the transmittance value at selected THz frequencies.
The results have been reported in Figure 9a-e, at frequencies 0.2, 0.7, 0.8, 1.0 and 1.3 THz, thanks to the absence of sample absorption spectral features in this spectral range (see Section S5 in Supplementary Materials).
The transmittance without gain (T real ) for each pixel in each image, can be obtained by using the following Equation (10): where T displayed is the transmittance for each pixel, displayed in Figure 9a-e, after applying a gain g. The value of g, for each image, was chosen in such a way that the maximum transmittance value becomes 100% transmittance. By increasing the frequency, the spatial features are better discriminated. g is reported under each Figure 9a-e. ~600 µm. Linen spectral images have been acquired. We collected transmittance images at different frequencies. Each pixel contains the transmittance value at selected THz frequencies. The results have been reported in Figures 9a-e, at frequencies 0.2, 0.7, 0.8, 1.0 and 1.3 THz, thanks to the absence of sample absorption spectral features in this spectral range (see Section S5 in Supplementary materials). The transmittance without gain ( ) for each pixel in each image, can be obtained by using the following Equation (10): The image corresponding to 0.2 THz is a uniform distribution of intensity pixel values. The low contrast and spatial resolution result in no distinguishable linen structure, due to its strong sub-diffraction nature. At 0.7 THz, a horizontal dark line starts to become visible. It is associated to the presence of a not resolved linen fiber (or more). Moving to high frequencies and increasing the spatial resolution, more linen fibers are identifiable, as it is fibers (black zones) are visible, while the white zones are the hole between linen fibers. Figure 9f shows the digital image and can be compared to the 1.3 THz image, that has the best resolution, among linen spectral images. The two vertical linen fibers are well resolved because spaced by big holes (>150 µm), while the three horizontal pairs are not resolved and only three dark regions are visible.
In order to obtain image de-blurring, we performed image deconvolution. With respect to Equation (2) As seen in the de-blurred Figure 10c, the dark regions, corresponding to the linen fibers, become better identifiable, thus an image quality enhancement is visible, resulting in a sharper image. In Figure 10d,e we report the intensity profile and its derivative, for both the original and deconvolved linen images. The red boxes in Figure 10a,c show the pixel lines selected for intensity profile and derivative plots. In Figure 10e D 1 and D 2 denote the differences between derivative maximum and minimum for original and deconvolved images, respectively. By computing the ratio R = D 2 /D 1 an estimation of the quality improvement can be given. The deconvoluted image shows a 1.6-fold quality improvement.
resolved and only three dark regions are visible.
In order to obtain image de-blurring, we performed image deconvolution. With respect to Equation (2)  As seen in the de-blurred Figure 10c, the dark regions, corresponding to the linen fibers, become better identifiable, thus an image quality enhancement is visible, resulting

Conclusions
In this paper we characterize a THz-TDS system based on PCAs. We evaluate the performances of the system such as its spectral range emission and spatial resolution. Moreover, we reconstruct the 1D and 3D profiles of our THz beam, with two different methods: hole and knife edge. Both methods gave results that are in great accordance. Finally, we present a very simple and cheap resolution test that can be done with a linen pattern sample. We demonstrated that our system has the capability to spatially resolve objects with dimensions from~1.3 to~0.3 mm in the frequency range 0.5-2.5 THz, with a resolution of~0.7 mm at 1 THz. Finally, with the help of the experimental PSF, we show that the THz image quality can be considerably improved.