## 4 Results  
### 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 \(\pm\) 6 \(\upmu\)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 \(\pm\) 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  
Figure 3: (**a**) Four _spectral_ images of the same 200 \(\upmu\)m pinhole for evaluating our imaging system point spread function (PSF). Figures (20 px \(\times\) 20 px, 3 mm \(\times\) 3 mm) are collected at four different frequencies: 0.6, 1.0, 1.6 and 2.5 THz. The scale bar and the spectral resolution in each image are 0.3 mm and 0.1 THz, respectively. (**b**) Full Width at Half Maximum (FWHM) estimation from spot intensity profile. Black: experimental data. Green: Gaussian fitting curve. Data shown for 1.0 THz image. Figure 2: Schematic visualization of Gaussian beam parameters: the beam waist in the focal region (\(W_{0}\)) and along the \(z\) propagation axis (\(W(z)\)), the Rayleigh distance (\(Z_{R}\)), the beam divergence angle (\(\theta\)). See text for details. 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)\)  
\[f(x)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x-y)^{2}}{2\sigma^{2}}} \tag{6}\]  
where \(\mu\) and \(\sigma\) are the Gaussian expected value and variance, respectively used as fit parameters. According to [47], the FWHM is:  
\[FWHM=2\sqrt{2 ln2}\sigma \tag{7}\]  
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)). \begin{table}
\begin{tabular}{c c} \hline \hline
**Frequency (THz)** & **FWHM (mm)** \\ \hline
0.6 & 1.06 \(\pm\) 0.02 \\
1.0 & 0.67 \(\pm\) 0.01 \\
1.6 & 0.42 \(\pm\) 0.01 \\
2.5 & 0.29 \(\pm\) 0.02 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: FWHMs calculated by the Gaussian fit, at different frequencies, following Equation (7). Figure 4: Spatial resolution vs frequency. In black: Full Width at Half Maximum (FWHM) calculated according to Equations (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.