Increasing the resolution of the reconstructed image in terahertz pulse time-domain holography

In this paper, we present a novel numerical approach for increasing the resolution of retrieved images of objects after their diffraction patterns are recorded via terahertz pulse time-domain holography (THz PTDH). THz PTDH allows for spectrally resolved imaging with high spatial resolution and does not require the fine alignment of complex optics in the THz path. The proposed data post-processing method opens up the possibility to reconstruct holograms recorded with spatially restricted THz detectors, and overcome the diffraction limit even for the lower-frequency spectral components. The method involves an iterative procedure of backward-forward wavefront propagation to simulate the field distribution beyond the initially recorded hologram area. We show significant improvement in both the object reconstruction and contrast across the whole spectrum, with qualitative resolution enhancement at lower frequency spectral components.

From traditional holography in the visible range, it is well known that an object can be retrieved from an arbitrary small hologram, with resolution inversely proportional to its size 48 . In lensless digital holography, an additional restriction is set by the number of pixels in the digitized hologram. The hologram size and number of pixels are even more vital in THz holography. On one hand, the sensor area cannot be immediately extended to cover a wider spatial angle, while on the other hand, the pixel size cannot be reduced since it is already equal or less than THz wavelengths, reaching from ∼100 µm to ∼3 mm.
Recently, for visible range holography, Latyshevskaia et al. proposed a self-extrapolation technique 49 which allows for overcoming the Abbe criterion and resulting resolution enhancement of the reconstructed objects by padding the hologram and iteratively filling the padded area with the numerically propagated wavefront. Theoretical justification (more detailed in 50 ) of the complex wave spatial extrapolation was discussed by J.L. Harris 51 , who has shown the uniqueness of the problem solution for given spatial spectrum distribution of the analytical signal on the basis of Whittaker Watson law "the spectrum of a size-limited object is always analytical" 52 and Guellemin law "a function of a complex variable is determined throughout the entire Z-plane from a knowledge of its properties within an arbitrarily small region of analyticity" 53 Later, this method was extended to phase retrieval 50,54 and inline narrowband THz holography 33 .
In this paper, we propose an adaptation of the above mentioned technique for resolution enhancement of broadband THz pulse time-domain holograms. This method is in high demand, due to the absence of large (in the numbers of wavelength) area THz detectors. We show that the enhancement of the holograms occurs only at wavelengths that yield at least one full interference fringe within the recorded hologram. To the best of our knowledge, this is the first demonstration of the method for broadband radiation.

Terahertz Pulse Time-Domain Holography
To record THz pulse time domain hologram one must first record the amplitude and phase of a collimated wavefront at the detection plane, then place the object of interest to record the diffraction pattern with the same detector. As a matter of fact it is a broadband inline holography with coherent detection provided by pulsed THz that gives information on both amplitude and phase distribution upon measurement.
There are two main techniques for coherent THz detection: (i) free-space electro-optical sampling (EOS) 55 or; (ii) detection in photoconductive antennas (PCAs) 1 . Each technique has its own advantages and disadvantages 56 , but they can both be implemented in synchronous and asynchronous 57,58 layouts, providing reasonably comparable results, the latter being more fast but at the same time requiring significantly more expensive setup. The PTDH resolution enhancement approach we offer can be applied to any coherent wavefront THz detector, and, as there are currently no array PCA detectors, we propose the following layouts for full wavefront EOS using either a synchronous, or asynchronous approach (Fig. 1).
The setup on the left is based on a classical EOS THz detection layout. The more advanced setup on the right uses a retroreflected probe IR beam, thus allowing for closer object-to-hologram plane distance, and therefore, a larger spatial angle that positively affects the resolution. Although the full wavefront is detected in both layouts, the size of the hologram is constrained by size of both the matrix detector and the electro-optic crystal. The detection sensitivity also imposes a limit on the hologram size, since commercially-available coherent broadband THz transmitters are not powerful enough to be effectively detected in larger cross-sections.
Consider the object  ν x y ( , , ) placed in the THz wavefront. After diffraction on the object, the incoming THz field E x y t ( , , ) is distorted into E x y t ( , , )  . Further propagation to the registration plane at the distance l forms the diffraction pattern ′   E x y t ( , , )  . By applying a Fourier transform, this pulse wavefront can be converted into the frequency domain ( , , ) is the Fourier transform of the hologram recorded in the time domain at the distance L from the object plane, and can be considered as a union of monochromatic holograms calculated for every frequency v.
Generally, we can consider the set of independent monochromatic objects  ν x y ( , , ) spatially digitized into M × M pixels to lead to a transverse size of D = (MΔx, MΔy), corresponding to pixel size of Δx × Δy. This set of objects produces a set of independent holograms ν   x y ( , , )  , each being N × N pixels and the overall transverse The number of both sub-objects and sub-holograms in the sets is equal to the number of frequency components in the spectrum 42 . Likewise, both object and hologram can be represented as sets of spatial-temporal field distributions, or spatial-spectral intensity charts. This concept is depicted in Fig. 2.
The THz hologram is a three-dimensional array of data recorded after diffraction on an object. This array carries information about the spatial-temporal, or spatial-spectral distribution of the THz field.
To reconstruct an object from a hologram, the recorded complex field is back-propagated into the object plane: 1 is the inverse wavefront propagation operator, here,  and  are the sets of all holograms and objects, correspondingly. Depending on the distance between the object and the hologram, it can be either a Rayleigh-Sommerfeld integral, 2D Fourier or 2D Fresnel transforms, or angular spectrum calculation 39,59 .

Iterative self-extrapolation algorithm for PTDH
The principle for object self-healing in the case of broadband THz PTDH is quite similar to the one described for monochromatic holography 49 . Proposed technique is based on the assumption that the diffraction pattern, recorded in the hologram, is created by waves that exist not only within, but also outside the registration area. Therefore, self-extrapolation of the digital hologram to the area outside of the experimentally recorded diffraction pattern can be performed. The principal scheme for the algorithm in the case of broadband THz PTDH is shown in Fig. 3. In general words, the hologram is back-propagated to reconstruct an object from the data 'as is' , as described in equation (2). Then, after the application of the mask that restricts the energy distribution in the object plane and helps to avoid the ghost image appearance, the reconstructed object wavefront is numerically propagated to form a synthesized hologram is the forward propagation operator, that, similarly to inverse operator −1  is chosen depending on the distance between object and hologram. Thus, we can introduce an iteration operator , describing the reconstruction of an object ∼ , and consecutive generation of a synthesized hologram ′ ∼  as follows: 1 When this operator is applied, one should keep in mind that any region of the object plane that is involved in calculations, but does not contain meaningful information, must be zeroed to avoid ghost images 60    After a reasonable number of iterations, the reconstructed object is recovered, along with its spectral features.

Results and Discussion
Since the spectrum of THz pulsed radiation is more that 1 octave wide, depending on the object type, the algorithm can be implemented in a different manner. For spectrally uniform objects, like binary amplitude opaque (e.g. metal) objects, reconstructed data obtained with one of the wavelengths can be used as a synthetic object in the next iteration for forward propagation simulation with a different wavelength. Thus, an enhancement across the whole spectrum can be achieved. On the other hand, for an object with numerous spectral features, each feature can be reconstructed independently, provided that the numerical aperture of the hologram gives the sufficient amount of wavefront data for reconstruction. Indeed, if the full contrast interference fringe is not recorded, the object will not be retrieved 49 .

Image reconstruction by low frequency components of the THz range.
Let us consider the problem of increasing the resolution of a binary amplitude object. In accordance with the Abbe criterion, the resolution of the reconstructed image for each spectral component is limited by the numerical aperture NA of the hologram, hence it depends on the distance between the object and the registration plane, as well as the linear dimensions of the hologram as Abbe where R Abbe is the resolution in m, c is the speed of light, D and l are the transverse size and distance to the hologram. Thus, after a size restricted hologram  ν   x y ( , , ) 0 is recorded, application of the iterative algorithm and hologram self extrapolation obtains the field ν ∼   x y ( , , ) i  , which has a larger transverse size, thus enhancing the reconstructed image resolution. Figure 4 illustrates the results of the algorithm applied to the binary object of the two apertures separated by ∼15 mm. The object itself, 25 × 25 mm 2 in size, is shown in the inset of Fig. 4(a). The central part of Fig. 4(a) shows the simulated 25 × 25 mm 2 hologram at 0.4 THz, recorded at 600 mm from the object. This frequency is selected, as in this frequency range the algorithm application demonstrates the highest effect. Indeed, if we calculate the initial resolution at 0.4 THz, it is equal to = R 18 bbe 0 A mm, and the result of the object reconstruction from this hologram is shown in Fig. 4(b). Then, this reconstruction is used to simulate the hologram , as it was described in Fig. 3. The wavefront in this hologram is wider than the area of the initially recorded hologram, hence the NA and the reconstruction resolution is higher ( = . R 4 5 bbe 1 A mm). The reconstruction after first iterative procedure shown in Fig. 4(d), already reveals the start of the two dots separation. After the 10th algorithm iteration, the theoretical resolution reaches = R 3 bbe 10 A mm, and the two dots in the reconstruction are clearly resolved. Holograms and reconstructed objects at further iteration steps (5th and 10th), are shown in Fig. 4(e,g) and (f,h), respectively.
Spectra of the reconstructed images cross-sections ( Fig. 4(i-l)) demonstrate the effectiveness of the proposed method. While slightly enhancing the reconstructed images in the higher frequency range, the crucial change is revealed in the spectral area between 0.35-0.45 THz (670-850 μm), where the wavelength is approximately equal to the detector pixel size (780 μm) and the theoretical resolution changes from mm upon iterations, i.e. from below to above the dots separation, and the dots become resolvable. At higher frequencies, the dots are being resolved even after the direct hologram reconstruction, and at longer wavelengths, the NA of the hologram is still not enough to cover at least one interference fringe and thus provide enough contrast in the object reconstruction. As a result, the apertures, although they are larger than the wavelength, can not be recovered through the proposed algorithm. Such a fast convergence of the algorithm, as compared with results in 49 , occurs due to longer THz wavelengths involved (hence smaller objects, holograms, and distances between them in wavelength units) and higher impact of every iteration step. Note that reconstruction of the image at the frequency of 0.42 THz require only 10 iteration, while at lower frequencies (e.g. 0.38 THz) more iterations should be performed for significant improvement of reconstructed image spatial resolution.
Application of the technique to spectrally selective object reconstruction. Next, we study a spectrally selective spatially inhomogeneous object (Fig. 5). As an object, we used opaque mask 44.8 × 44.8 mm 2 in size with 'P' , 'T, 'D' and 'H' letters transparent in different spectral windows.
The spectral-selective amplitude transmission of an object is given by the function: Here R(x, y), G(x, y) and B(x, y) are spatial normalization coefficients of the amplitude for the Gaussian distribution; ν ν ν , , The transmission spectra are shown in Fig. 5(c). The letters P, T and D in Fig. 5(a) are transparent in the lower (red), medium (green) and high (blue) frequency ranges, respectively. Meanwhile, the letter H is transparent in all three frequency ranges. The size of the square object in Fig. 5(a) Table 1 present the results of the iterative object reconstruction. Both spectral and spatial features eventually get resolved after ∼10 iterations, and the object gains larger area of fine reconstruction, which initially was restricted by the smaller size of the hologram and very low scattering introduced by the objects. As a result, the mean root square error (MRSE) of reconstruction quality (normalized error between reconstructed image compared to initial object) drops from 0.66 to 0.34 (see Table 1). Due to the fact that MRSE depends mostly on the general shape of the object rather than on fine structures of the object (in this case sharp edges of the letters), significant decrease of MRSE can be observed right after the first iteration, which provides overall shape of the letters. On the other hand, addition of higher spatial frequency components during further iterations, which results in increase of letter edge reconstruction quality, does not introduce significant contribution to decrease of MRSE.

Conclusion
In this paper, we proposed the solution for resolution and field of view enhancement in THz pulse time domain holography (THz PTDH). THz PTDH allows for spectrally resolved images to be obtained without any additional optics in the THz beam, as the signal is coherently detected after diffraction on the object of study. The resolution enhancement method involves an iterative procedure of backward-forward wavefront propagation, with numerical extension of the hologram area. The method is especially important for THz holography, as the detectors are usually restricted in size, and the relatively large wavelength means that a very low number of interference fringes can be recorded. We have shown the method to enhance both the resolution and the field of view for the reconstructed images after performing only 10 iterations. We successfully applied the method for both spectrally uniform and spectrally selective objects, and revealed significant enhancement at all frequencies that produce at least one interference fringe in the diffraction pattern of the recorded hologram. At the frequencies between 0.35-0.45 THz, the Abbe criterion for the object and hologram under study is reduced from the value of 23 mm (above the distance between the apertures used as an object) at the first iteration, to the value of 3 mm (well below the distance between the apertures) at the last iteration. Spectrally selective test object also revealed significant reconstruction quality enhancement resulting in double reduction of reconstruction mean root square error. Thus, the proposed method can be can be successfully applied to all kinds of pulsed THz holographic imaging tasks, a well as to spatio-temporal metrology of broadband complex wavefronts 59,61 . Field of view broadening, in all these tasks, when THz detection is performed with relatively small crystal, will push the boundaries of the method experimental applicability.

Methods
Numerical simulations were performed with in-house software, written in National Instruments LabView development framework. The software follows mathematical procedures described in this paper, and our preceeding works 39,59,61 .