Single-scan multiplane phase retrieval with a radiation of terahertz quantum cascade laser

Terahertz phase retrieval from a set of axially separated diffractive intensity distributions is a promising single-beam computational imaging technique that ensures the obtention of high spatial resolutions and phase wavefronts, but remains restricted by time-consuming data acquisition processes. In this work, we have adopted an approach, relying on the radiation of a quantum cascade laser and the implementation of an express single-scan measurement of intensity distributions through the continuous on-the-go displacement of a high-sensitivity antenna-coupled microbolometer sensor array. In addition to the simplicity of this practical implementation and the minimization of measurement times, such an approach overcomes the problem of preliminary optimal selections of transverse intensity distributions used in the iterative phase retrieval algorithm and guarantees the required data diversity for high-quality wavefront reconstruction.

A multiple-plane phase retrieval is an efficient approach, which relies on the measurements of a set of axially separated transverse intensity distributions. Similar to the competing techniques, it also makes use of processing algorithms adapted from methodologies developed in the visible frequency range [42,43]. The most popular algorithms are single-beam multiple intensity reconstruction (SBMIR) [53] and the transport of intensity equation [48]. This multipleplane approach is characterized by the extreme simplicity of the object illumination setup and guarantees better convergence compared to single-plane phase retrieval algorithms, but the recording of the diffraction patterns often remains time-consuming.
Another problem that one has to face in multiple-plane phase retrieval in the THz frequency range is the existence of a narrow range of recording distances to the object for the optimal positioning of the sensor, from which it is possible to ensure the highest quality of the reconstructed images [43]. These optimal positions depend on the features of the diffraction process and can be generalized in terms of sets of dimensionless Fresnel numbers. It has been shown that registering intensity distributions with a too-small longitudinal pitch does not guarantee sufficient diversity in the data structure required for the convergence of the iterative algorithm [43]. Oppositely, the registration of transverse intensity distributions with excessively large longitudinal intervals does not lead to sufficiently high resolutions on the reconstructed images since the content of high spatial frequencies fades away for far-field diffraction. In general, the optimal longitudinal step when registering diffraction patterns separated along the optical axis depends on the diffraction zone, which is conveniently characterized by the distance between the object and the nearest intensity measurement plane.
In addition, it should be noted that the THz frequency range differs greatly from the visible wavelength range through the loss of the paraxial nature of the diffraction processes, which is usually caused by the small transverse dimension of a typical THz beam, with respect to the operating wavelength [20]. Thus, the choice of the optimal number and arrangement of measurements planes presents certain complexity and often have to be empirically assessed, and carried out through several series of experimental measurements or preliminary numerical simulations, hence remaining a time-consuming anticipation process.
And finally, success in solving practical problems of THz PI is guaranteed by the choice of a suitable combination of sources and detectors of radiation in this range of the electromagnetic spectrum. Early work on THz PI [4,[23][24][25]28] used single-pixel detectors in the raster scan mode of operation due to the lack of matrix detectors. As array detectors became available, a larger variety of THz PI works initiated such implementations [54][55][56][57][58]. Nevertheless, the limited sensitivity of the THz radiation matrix detectors, along with the recurrent mismatches between the immutable pixel pitches and wavelengths of the detected THz radiation, do not allow the adaptation of existing solutions for an arbitrary combination of source and detector. To illustrate this point, a previous experiment demonstrated a multiple-plane phase retrieval setup in a reflection configuration based on a sequential raster scan using a single-pixel Shottky diode sensor, aiming to assess the diffraction pattern generated from up-converted Gunn diode at the wavelength of ∼ 1 mm [53]. Other works, contrarily, made use of matrix detectors for the assessment of the set of intensity distributions in the diffraction field [48,59], coupled to a high-power optically pumped far-infrared gas laser as the radiation source.
Thus, the given work aims to cover such technological shortcomings for multiple-plane phase retrieval, namely by overcoming (1) the currently limited demonstrations of source-detector combinations, the (2) time-consuming data acquisition, and the (3) necessity of preliminary optimizations of the experimental conditions. More specifically, we successfully demonstrate single beam multiple-plane iterative phase retrieval through the radiation of THz quantum cascade laser (QCL) and a real-time recording via a highsensitivity antenna-coupled microbolometers camera. This experimental setup allows the express assessment of all the intensity distributions necessary for successful wavefront reconstruction in a unidirectional scan, with continuous acquisition implemented during a single longitudinal displacement of the array sensor. Thanks to the sensitive detection hardware, the implemented method is convenient and ensures an efficient way for data acquisition in on-the-go (OTG) mode. Hence, it eradicated the necessity of a preliminary choice concerning the optimal number and arrangement of measurements planes. It then guarantees the required data diversity for high-quality wavefront reconstruction and opens up the possibility for further automatic adjustment for the parameters of the iterative data processing.

Experimental implementation
In lensless computational diffractive THz PI, the bare bolometric array camera ('C' in Fig. 1), deprived of the imaging lens, can be employed for the simultaneous two-dimensional intensity distribution measurement in the transverse planes of the diffraction field. In a straightforward transmission configuration, such a recording unit should be complemented with the illumination source, required beam-shaping optics, sample holder (denoted as 'O' on Fig. 1), and the motorized translation stage ('TS') for the measurements of the diffracted intensity profiles, I z k (x, y) , separated along the optical axis. In this work, Lytid's TeraCascade 1000 QCL source, with an emission power of 1.3 mW at 2.5 THz, has been employed in combination with its integrated autoalignment unit ('AAU' in Fig. 1) to provide a collimated illumination beam through 1" optics [2,60].
Based on beam shaping optics and deflecting mirrors, Lytid's auto-alignment unit performs an automated source selection on the QCL multi-chip (there can be up to 6 chips in the Teracascade source). Then, from the diverging QCL's emission profile, it ensures the provision of a collimated beam, with a wavefront profile close to plane, to obtain a highly homogeneous illuminating wave phase, especially suitable for such applications. The example of the intensity distribution of the collimated beam captured without the object is shown in Fig. 1. Such a limited illumination diameter finds itself suitable when ranging up to 2.5 THz, with respect to lower frequencies implementations [53] that require larger inspection fields. Indeed, an obvious setup and sample rescale should be expected relative to the wavelength of operation.
In this transmission configuration of experimental layout, a polypropylene (PP) phase object (see Fig. 2) with thickness of 1.59 mm in unpatterned areas and a sample refractive index of n PP = 1.5151 at 2.5197 THz [41], has been investigated. Its specific topographic features display details of the order of the expected lateral resolution. Two regions of interest (RoI) of the sample were studied in this work: plastic identification code ('PIC') and 'UP' with grooved and extruded patterns, respectively, as shown in detail in the insets of Fig. 2.
A preliminary experiment has been performed using a first-generation INO IRXCAM camera ( 288 × 384 pixels matrix with 35 μ m pitch), which displayed stability issues over the recording operations, hence limiting the measurement repeatability. Additionally, a drastic lack of sensitivity, coupled with the milliWatt illumination level, constrained the intensity recording to diffracted low spatial frequencies.
The subsequent SBMIR wavefront extraction then only led to relevant reconstructions on restricted portions of the illumination field. The tenfold increase in sensitivity with the use of an I2S camera (240 × 320 pixels matrix with 50 μ m pitch and detection rate 25 frames per second (FPS)) has then been employed for the suited recording of the intensity diffraction field data sets.
A prior stop-motion (SM) capture of the diffraction field, I z k (x, y) , was initially approached. The SM acquisition process is the conventional intensity distributions measurement mode, which consists of the sequential movement of the sensor for certain longitudinal steps Δz and the data recording in between these movements at the moments when the detector is stationary. It became widespread because it ensures the perfect mechanical stability of the sensor when recording, while potentially allowing for multi-frame averaging to improve the measurement SNR on each captured diffraction pattern (further, such an approach will be referred to as SM-A). As a rule of thumb, such capturing procedure takes to several minutes and the exact acquisition time depends on the desired number of frames used for averaging. In the case considered in this paper, we used a longitudinal pitch Δz = 2 mm, as featured in Fig. 3a, c and at each position, a set of roughly 35 images of the given diffraction patterns were captured for subsequent manual averaging. Due to the non-automated procedure, the time spent between every measurement for the sensor displacement varied. On average, the detector movement time interval was 2.184 ± 0.08 s. Therefore, the acquisition speed in the SM-A mode can be estimated as the number of pixels in the resulting averaged image of a certain measurement plane, with respect to the total movement and capturing time for the whole series of images being averaged. The result of the calculation is presented in Table 1. When not considering any additional averaging, the time spent on registering data in the SM Fig. 1 Diagram of the experimental implementation for real-time phase retrieval measurement in transmission configuration (a) and the example of a formed far-field beam pattern, captured without object (b) Fig. 2 Depiction of the PP Phase sample, investigated in this work, featuring two topographic RoI mode theoretically can be neglected. We will refer to this option as SM-S, and the estimated acquisition speed for such approaches are presented in Table 1, notwithstanding we did not implement it in our experiments.
Featuring a continuous displacement of the sensor, the OTG data acquisition was also implemented. In this operating method, the movement of a motorized translation stage was initiated along with synchronization of the measurement process. Data acquisition with a continuous displacement of the sensor at the speed of 10 mm/s was defined as an optimal displacement speed in this implementation. It requires a recording time below 1 second for a longitudinal displacement of a few millimeters. Accounting for the recording rate of 25 FPS, a Δz = 400 μ m longitudinal recording pitch is achieved for a raw recording data-set, as depicted in Fig. 3b, d. The acquisition speed for OTG mode then reaches 1.92 × 10 6 pix/s. Frame decimation in the diffraction profiles stack for the selection of an optimized input data-set for the SBMIR procedure can then be performed and refined as a post-processing operation.
The above-mentioned estimations of acquisition speed are then to be compared with the acquisition speed for the empirical single-pixel recording by raster scanning, implemented in previous work [53] (RS in Table 1).

Methods
From such diffraction fields recorded data-sets, displayed in Fig. 3, a square 240×240 pixel selection of the collected set, centered on the sample portion of the interest as well as the illumination beam, is cropped as to fit the subsequent algorithmic approach.
Namely, the straightforward sorted iterative SBMIR procedure features a strained propagation between measurement planes. The subsequent wavefront distributions are obtained via the employment of the free-space propagation kernel, H(f x , f y ) , offered through the Fourier optics mechanism and employed as follows:  where z is the arbitrary distance between the two subsequent planes, U(x, y, 0) and U(x, y, z) the complex amplitude of the wave front at the source plane and observation plane respectively. F and F −1 represents the direct and inverse Fourier transform operators respectively. Regardless of any approximation, the general transfer function for free-space propagation over a distance z, H(f x , f y ) , can be expressed as: where f x = k x 2 and f y = k y 2 are the spatial frequency components in the Fourier space, along the � ⃗ x and � ⃗ y directions respectively, is the wavelength.
For each considered measurement plane, a constraint in amplitude is then applied as to fit the intensity diffracted field measurement such that |U z k (x, y, where U z k is the estimated complex field. This operation will then lead to the convergence of the wavefront phase profile along this iterative process. Following this approach, a field focalization in the object plane can then be performed as to retrieve the object profile. Nevertheless, in such a configuration, a precise mechanical assessment of the first measurement plane distance, z 0 , is impracticable due to the obstruction of the camera mechanical structure, as well as the unreachable sensor plane, lodged U(x, y, 0) , behind a protective window. Such a direct focused reconstruction on the object plane, following each SBMIR step, remains then inconceivable. The wavefront reconstruction has then to be performed at a raw estimated object plane position and further focused in a subsequent effort.

Results and discussion
The focused coherent field reconstructions, emerging from the SM-A recorded data-set and the proposed optimized OTG recording are displayed in Figs. 4 and 5 for two object's RoIs: PIC and UP respectively. Figure 4 relates to the inspection of the PIC RoI through different SBMIR and data collection approaches. Namely, amplitude and phase images shown on Fig. 4a, d, correspondingly, were recovered from the first m = 5 intensity distributions of the dataset, obtained in SM-A mode, with the following parameters: z 0 = 11.9 mm, Δz = 2 mm. The use of only the first five input diffraction intensity distributions ensured the optimal quality of the reconstructed images while the inclusion of additional intensity distributions from more distant planes led to a worsening of the reconstruction quality, due to the lack of high spatial frequencies in that region. In addition, the single-frame captured for each of the measurement planes was extracted from the raw SM-A dataset (that is images before the averaging procedure) to simulate data recorded in SM-S mode. Negligible improvements have been witnessed when averaging over multiple still frames due to the adequacy of the emission power level with respect to the camera noise equivalent power. Fig. 4 Resulting numerical wavefront reconstruction, for amplitude (top rows) and phase (bottom rows) profiles, for 'PIC' RoI after 30 SBMIR iterations. The left column includes the data obtained using SM-A recording with Δz = 2 mm, while the middle and right columns, respectively, display wavefronts recovered using OTG recording with respective longitudinal pitches of Δz = 2 mm and at the optimal value of Δz = 2.4 mm In the latest proposed OTG configuration, the 25 FPS detection rate of the camera allows for a continuous displacement of the planar array, while preventing any blurring effect during the sensor motion in the input data sets. Indeed, high requirements to the mechanical stability are critical in phase retrieval, implemented with the radiation of visible frequency range where the wavelength is less than the possible detector deviations caused by vibrations. In the THz frequency range, the radiation wavelength is much greater than the detector possible offsets range. That is why data recorded by a sensitive detector in OTG mode do not provide any deterioration of the reconstructed images. Herewith, at the previously mentioned frame rate, the achievable longitudinal pitch, Δz = 400 μ m, remains beyond the required optimal value of recording spacing. It validated again, that, considering the smallest longitudinal pitch in the collected dataset does not necessarily lead to an optimal quality of the reconstructed images. This aspect has been extensively discussed in [43] as algorithmic optimizations have been assessed.
From the extensive OTG recordings, sub-datasets fitting the previous parameters, namely Δz = 2 mm, have been employed to reconstruct the wavefront, with amplitude and phase characteristics displayed in Fig. 4b, e. With respect to the slower SM-A recording, it stands out that the suitability of the proposed OTG acquisition mode, complemented with the collection of a decimated dataset selection, leads to the reconstruction of a slightly sharper coherent phase profile while a similar lateral resolution is obtained.
The versatility granted through the extensive data-set OTG capture is further highlighted in Fig. 4c, f. Indeed, further investigations demonstrated an optimal numerical longitudinal pitch of Δz = 2.4 mm, allowing to reach a representative wavefront recovery with once again z 0 = 11.9 mm and m = 5 . Through the provision of such an exhaustive initial OTG input data, the optimal parameter set can then be empirically found in a post-processing task, preventing repetitive data registration to reach peak performances.
To endorse once gain the suitability of the proposed OTG approach, Fig. 5 confronts several processing configurations on the resulting reconstructions for the 'UP' RoI. More specifically, contrary to the previous configuration, a further investigation has been done, where the first intensity distribution was selected at distance z 0 = 21.2 mm. Such a more distant configuration has been motivated by the presence of practical cases where the location of the detector array close to the sample remains impracticable, for example, due to mechanical obstructions [61]. Similar restrictions take place in the case of the phase retrieval setups, operating in the reflection mode [53].
Since in the second experiment intensity, distributions were registered further from the object at the different diffraction subband, to reconstruct an image with the highest achievable quality, the optimal longitudinal pitch also should be reconsidered. Specifically, Fig. 5a, d and b, e display the reconstruction improvement when ranging from Δz = 2 mm to Δz = 4 mm through the initial SM-A capturing scheme, when considering m = 5 recording planes. However, due to the poor longitudinal resolution provided through this tedious and time-consuming acquisition process, no finetuning of this parameter of interest can be expected without additional experimental measurements. The OTG approach nevertheless provides a suitable alternative to such limitations. Indeed, an optimal wavefront, depicted in Fig. 5c, f,   Fig. 5 Resulting numerical wavefront reconstruction, for amplitude (top rows) and phase (bottom rows) profiles, for 'UP' RoI after 50 iterations. The left and center columns include the data obtained using the SM-A recording method with Δz = 2 mm and Δz = 4 mm respectively while right column displays wavefronts recovered using the OTG recording method with an optimal longitudinal pitch of Δz = 3.2 mm has been reconstructed from a decimated input data-set featuring Δz = 3.2 mm.
Thus, in the light of the discussed results, it can be concluded, that the short recording longitudinal pitch, obtained through the high acquisition speed is especially suited as a prospective solution to be put to good use for the selection of an optimal SBMIR input data set, while allowing for fast full diffraction field assessment, under 1 second. Here, among this full recording, a numerical decimated set with Δz = 2.4 mm longitudinal pitch has been found to be optimal for the 'PIC' sample, hence considering 1 out of 5 400 μ m discreet diffraction profiles. The further assessment of the diffraction field led to an optimal longitudinal pitch of Δz = 3.2 mm for the 'UP' RoI.
Nevertheless, independently of the recording approach, those results display wavefront inhomogeneities in their phase and amplitude profiles. Those features can be devolved to the probe beam wavefront inhomogeneities, even though the careful beam shaping of the auto-alignment unit ensures the generation of a collimated plane wavefront from the diverging beam profile of the emission QCL. Indeed, the beam profile differs from Gaussian one, as shown in Fig. 1. These inhomogeneities can include wavefront singularities, negatively affecting the construction of a surface topography map. Nevertheless, knowing the spatial coordinates of the area with the most pronounced parasitic wavefront curvature, it is sometimes possible to implement a numerical correction procedure (Fig. 6a). However, it is not always convenient. Therefore, at some phase maps the impact of the parasitic wavefront curvature may remain uncorrected, leading to errors in a surface profile (Fig. 6b). Topographic variations between the unpatterned and grooved areas for the 'PIC' RoI and between the background and the extruded letters 'UP' are about 200 μ m and 100 μ m, respectively.
Beyond the addressed coherent imaging perspectives, the employment of such real-time recording techniques to provide coherent wavefront distribution can as well be put to good use for source characterizations and far-field wavefront sensing. It would indeed overcome the limited resolution of the sub-samples Hartman method [31][32][33][34], tied to the hole-punched mask and sensor pitch. In the case of QCL radiation, the wavefront of which still contains singular points after collimation, multiple-plane phase retrieval the can be easily implemented. However, in the case of the smooth, non-diffracting wavefront, the basic SBMIR algorithm became inefficient for wavefront sensing. To deal with this problem a small controlled disturbance must be introduced into the examined wavefront, which can be removed afterward from retrieved phase distribution, as was demonstrated in [53].

Conclusion
The tedious multi-plane data collection process, required for the SBMIR algorithmic approach remains a major drawback of such iterative phase retrieval procedures. Through a simplified transmission configuration, featuring a high-power QCL adequately coupled to a highly sensitive camera sensor, a real-time recording with the continuous displacement of the sensor allows OTG acquisition of the full diffraction field, eligible for high-quality reconstruction. Such a fast recording process, ensures the collection of an exhaustive dataset that provides substantial flexibility in the postprocessing extraction steps of the wavefront. Namely, the importance of the selection of the numerical longitudinal pitch, Δ z , on the reconstruction quality has already been highlighted and the OTG recording allows this parameter to undergo an optimization in post-processing operations, preventing time-consuming iterative experimental tasks.
Beyond the scope of this work, a few limitations remain discernible on the wavefront reconstructions. Namely edge effects, induced by the limited sensor dimensions are impacting the wavefront uniformity. Through algorithmic extrapolation procedures [16,62], or experimentally, with the mapping of a wider imaging field through the lateral displacement of the camera [27,63] combined with the longitudinal multi-plane OTG recording, the enhanced retrieval of complex wavefront profiles would be expected.