# Computational terahertz phase imaging using a sequence of multi-plane intensity measurements

Vivek Kumar, Pitambar Mukherjee, Frédéric Fauquet, Kedar Khare, Sylvain Gigan, Patrick Mounaix

► To cite this version:

Vivek Kumar, Pitambar Mukherjee, Frédéric Fauquet, Kedar Khare, Sylvain Gigan, et al., Computational terahertz phase imaging using a sequence of multi-plane intensity measurements. APL Photonics, 2025, 10 (9), 10.1063/5.0283596. [hal-05291507](https://hal.science/hal-05291507v1)

HAL Id: hal-05291507

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

Submitted on 1 Oct 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
# Computational Terahertz Phase Imaging using a Sequence of Multi-Plane Intensity Measurements

Vivek Kumar<fnref n="1" />,<fnref n="2" />*, Pitambar Mukherjee<fnref n="2" />, Frédéric Fauquet<fnref n="2" />, Kedar Khare<fnref n="3" />, Sylvain Gigan<fnref n="1" />, Patrick Mounaix<fnref n="2" />

**Abstract:** Advancements in terahertz (THz) technology have substantially propelled the capabilities of imaging systems; however, retrieving the complex phase distribution of arbitrary test objects within the THz spectral range continues to face fundamental and practical challenges. In this context, interferometry-based imaging schemes have been widely adopted for THz phase imaging. However, the efficiency of such techniques is fundamentally constrained by their reliance on phase-stable reference beams, rendering them particularly susceptible to noise and inherent imaging system instabilities. In this work, we address such limitations and experimentally investigate two non-interferometric imaging configurations integrated with a computational phase retrieval framework, enabling accurate recovery of complex-valued wavefronts from intensity-only data. The first configuration employs a lens-less imaging approach based on recording axial diffraction patterns at multiple planes. The second configuration captures defocused intensities near the focal plane of a 4f imaging system. Our proposed computational framework formulates the phase retrieval task as a nonlinear least-squares error optimization problem. We employ an accelerated second-order optimization strategy by iteratively computing both the gradient and the Hessian of the cost function, which outperforms previously demonstrated multi-plane phase reconstruction approaches. Our imaging framework marks a significant step forward in computational THz imaging, with promising implications for diverse applications such as non-destructive testing for material characterization, biomedical diagnostics, and security screening.

**Keywords:** terahertz phase imaging, nonlinear convex optimization, multi-plane phase retrieval

**Introduction:**

Phase imaging in the THz spectral range has opened up a wide array of potential applications in biomedicine, security, and quality control, thanks to the non-destructive and non-ionizing nature of THz radiation<fnref n="1" />⁻⁵. Generally, the spatial phase information enhances imaging contrast for transparent samples, offers essential insights such as the relative optical depth profile, internal density variations, and encodes topographic details of the sample that are otherwise difficult to access⁶. However, despite these advantages, THz phase imaging encounters several fundamental and practical challenges⁷. In contrast to phase imaging in the optical frequency domain, where well-characterized illumination wavefronts and highly sensitive cameras enable phase retrieval, THz phase imaging is often constrained by factors such as low spatial coherence of illumination, lack of bright sources, limited dynamic range of image sensors, and the constrained space-bandwidth product of current state-of-the-art detectors. Consequently, considerable research efforts have been invested to develop robust THz phase imaging methodologies, particularly leveraging modern state-of-the-art time-domain spectroscopy (TDS), which provides not only amplitude and phase, but also ultrafast temporal and spectral information of the imaged object⁸⁻¹². Within this framework, several imaging modalities, such as near-field imaging¹³,¹⁴, Fourier synthetic aperture imaging¹⁵, ghost imaging¹⁶⁻¹⁹ and its application to imaging through complex media²⁰⁻²², have emerged, offering unique capabilities that are often unattainable in other spectral domains. Despite THz-TDS remaining the dominant technique for extracting spectral phase information, it inherently depends on

<footnote n="1">¹Laboratoire Kastler Brossel, ENS-Université PSL, CNRS, Sorbonne Université, Collège 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">³Optics and Photonics Centre, Indian Institute of Technology Delhi, New Delhi 110016, India.</footnote>

<footnote>*Corresponding author email: vivek.kumar@lkb.ens.fr</footnote>
femtosecond laser pulses, which impose considerable system complexity and cost and limit the achievable brightness of THz sources. Additionally, due to the absence of matrix detectors, the TDS technique requires raster scanning of samples, significantly extending data acquisition and processing times. Furthermore, maintaining stable THz sources typically requires stringent, controlled laboratory environments.

Recent technological advancements have fostered the development of commercial THz instruments with enhanced performance, improved reliability, and greater ease of operation. In this context, continuous-wave (CW) THz radiation can now be efficiently generated and detected using modern electronic devices that are compact, cost-effective, and straightforward to implement. However, most THz-CW imaging techniques are limited to capturing only the complex wavefront intensity distribution from an imaging sample, lacking information about the spatial phase distribution. In response to emerging challenges, CW-THz phase imaging has been extensively investigated, leading to the implementation of holographic imaging modalities both in off-axis²³⁻²⁶ and in-line²⁷⁻³⁰ configurations. Since the initial demonstrations of digital holography employing 2D THz image sensors, THz digital holography has advanced rapidly, benefiting from both hardware innovations and algorithm enhancements³¹,³². For instance, techniques such as dual-frequency reconstruction have been used to address 2π phase ambiguity³³, multiple plane imaging has been adopted to suppress image artifacts³⁴,³⁵, lateral shifting of the image sensor array has been employed to create a larger detection area²⁴,³⁶, sub-pixel shifting has enabled pixel super-resolution²⁸,³⁷, and physics informed deep learning for phase retrieval³⁸. Despite these remarkable developments in holography-based imaging techniques, challenges remain, particularly regarding reference arm phase stability, which can introduce undesirable phase noise and negatively impact image quality. Moreover, THz digital holography faces additional constraints when imaging multi-layered or highly lossy objects²⁶. Along the same direction, the Ptychography imaging approach has been adopted in the THz spectral range³⁹,⁴⁰. While THz Ptychography provides a reference-less and lens-less solution to the phase problem, it is often hindered by prolonged data acquisition times and slow algorithm convergence⁴¹. Quite interestingly, THz diffractive imaging employing multi-plane phase retrieval algorithms has been actively explored as a compact solution to the phase problem⁴²⁻⁴⁵. This approach leverages the sequential measurements of diffraction patterns at uniformly spaced propagation distances and uses a numerically backward-forward wave propagation scheme to iteratively recover the spatial phase information of the test sample. However, the multi-plane phase retrieval technique is fundamentally constrained by the requirement for sufficient redundancy and diversity in the recorded intensity measurements. Achieving such redundancy typically necessitates a larger number of measurement planes, which, in turn, exacerbates the accumulation of noise across planes. As a result, the quality of image reconstruction degrades and limits the broader applicability of the multi-plane phase retrieval approach.

Typically, in stark contrast to interferometric imaging methods, reference-less phase imaging involves recovering complex-valued elements of the object field directly from intensity measurements. As it pervades a diverse array of applications, prompting many researchers to focus on the development of phase-retrieval algorithms⁴⁶⁻⁴⁹. Classical approaches encompass a wide range of techniques, including generic gradient-descent routines, which become closely associated with the next generation of phase reconstruction algorithms and push the limit of conventional imaging⁵⁰. In this context, the gradient descent method approximates the search direction by utilizing the first-order derivative (gradient) of the cost function⁵¹. However, because the optimization is non-convex, gradient descent can easily become trapped in local minima or saddle points. Furthermore, its performance in terms of convergence rate can be limited, and it often fails to provide accurate phase reconstructions, particularly when dealing with low signal-to-noise ratio (SNR) intensity measurements⁵². To address such challenges, various strategies are available, ranging from conjugate gradient⁵³ and Nesterov acceleration⁵⁴,⁵⁵ to second-order methods⁵⁶. Notably, employing the second-order derivative (Hessian) of the cost function as a search direction improves the phase reconstruction fidelity⁵⁰,⁵⁶. Nevertheless, calculating the second-order derivative of the cost function can be computationally intensive in practice. Although several techniques have been proposed for implementing second-order optimization, they often prove impractical due to their significant memory requirements for storing and computing the Hessian or its inverse⁵⁰,⁵⁶,⁵⁷. Moreover, second-order optimization techniques have been explored in the context of phase retrieval within the visible spectrum,
however, such approaches typically operate under extended illumination conditions and pursue the joint estimation of the complex object field and the source distribution through alternating update schemes, analogous to a blind deconvolution framework<sup>58-60</sup>. Following this, an interesting question is whether the second-order optimization schemes can be exploited with a better strategy to develop efficient and noise-robust phase recovery algorithms in the THz spectral range.

In line with this vision, we introduce a non-interferometric THz phase retrieval technique integrated with a computationally assisted phase retrieval routine. We experimentally demonstrate two imaging configurations employing a CW THz quantum cascade laser (QCL) source operating at 3.5 THz and a direct intensity detection scheme using a microbolometer array detector. In the first configuration, we acquire multiple diffraction intensity patterns at different axial planes in a lens-less arrangement. In the second, we record various defocused intensity patterns near the focal plane of a 4f imaging system. In both the non-interferometric experimental approaches, we utilize recorded multiple intensity patterns transmitted through the sample to establish a nonlinear convex optimization framework built upon the angular spectrum formulation to describe the THz wave propagation model. Specifically, we use a second-order optimization method, utilizing the functional gradient and Hessian of the cost function to iteratively reconstruct the spatial phase information. To facilitate this, we derive analytical expressions to directly evaluate the function gradient and Hessian of the cost function, enabling its efficient computation through matrix-vector multiplication and fast Fourier transform operations. Consequently, our method circumvents the additional memory requirements typically associated with computation of functional gradient and Hessian of the cost function. As a result, our approach delivers a reliable, accurate, and rapid solution for addressing the second-order optimization problem, a key element in our phase imaging technique.

## Methodology

### THz phase retrieval using a second-order optimization framework

We present computational phase recovery of an arbitrary test sample using a second-order optimization approach. Figure 1 shows the schematics of experimental settings, where a collimated THz beam illuminates the sample placed at the object plane, and the transmitted two-dimensional complex wavefront intensities are captured at multiple axial planes using a THz micro-bolometric array camera mounted on a motorized translation stage. In this context, we collect intensity patterns in two imaging configurations. As shown in Figure 1a, we employ a lens-less computational diffractive THz phase recovery approach by measuring transmitted diffraction patterns directly at various planes along the optical axis. In the second imaging configuration, shown in Figure 1b, we collect the transmitted scattered wave field intensities from the sample near the focal plane of a telescopic arrangement of two lenses, forming a classical 4f imaging system. In both experimental schemes, our objective is to reconstruct the complex field of the imaging sample from a series of intensity images. In such a direction, we begin by introducing the forward model, which mathematically encapsulates the underlying physical process. To this end, we employ the angular spectrum formulation<sup>6,61,62</sup> to describe the THz wave propagation model, which, in its continuous form, is expressed as:

$$E(x, y, z_{\pm}) = \mathcal{F}^{-1} \left\{ \mathcal{F}\{E(x, y, z_0)\} \exp \left( ikz_{\pm} \sqrt{1 - (\lambda f_x)^2 - (\lambda f_y)^2} \right) \right\} \quad (1)$$

where $\lambda$ represents the wavelength of incident THz radiation, $k = 2\pi/\lambda$ is the associated wavenumber, and $z_{\pm}$ refers to the forward or backward propagation distance from the reference plane at $z_0$. The terms $E(x, y, z_0)$ and $E(x, y, z_{\pm})$ correspond to the complex-valued electromagnetic fields at the input and propagated planes, respectively. The variables $(x, y)$ and $(f_x, f_y)$ denote the spatial and spatial frequency coordinates, while $\mathcal{F}$ and $\mathcal{F}^{-1}$ represent the forward and inverse Fourier transforms, respectively. Following the forward model outlined in equation 1, a series of intensity measurements is obtained by axially translating the THz camera along the optical z-axis,

$$I_n(x, y, z_n) = |E_n(x, y, z_n)|^2 \quad \forall n = 1, 2, 3, \dots, N \quad (2)$$
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Second-order optimization-based THz phase imaging. (a) Schematic of THz lens-less phase imaging using multi-plane intensity measurements. The experimental configuration illustrating the acquisition of multiple diffraction patterns via axial translation of the THz camera, followed by phase reconstruction using a second-order optimization routine. (b) THz phase retrieval from a stack of defocusing intensity measurements. Schematic of the experimental settings showing the collection of defocusing intensity patterns near the focal plane within a 4f imaging configuration via axial translation of the THz camera, followed by phase reconstruction using a second-order optimization framework.</figcaption></figure>

where $n = 1, 2, ... N$ are the number of intensity patterns and $I_n(x, y, z_n)$ denotes the $n^{th}$ recorded intensity distribution at propagation distance $z_n$. Further, following the standard approach, we discretize the forward model using the matrix-vector formalism common in computational optics. In this framework, each intensity image $I_n(x, y, z_n)$ is vectorized into a column vector $\mathbf{I}_n \in \mathbb{R}^R$, and the object field $U(x, y, z_0)$ is similarly flattened into a complex-valued vector $\mathbf{U} \in \mathbb{C}^R$, where $R$ are independent spatial segments or total number of pixels in the image. The resulting discrete forward model takes the form:

$$ \mathbf{I}_n = |\mathcal{D}^\dagger \Psi_n \mathcal{D} \mathbf{U}|^2 \quad (3) $$

Here, $\mathcal{D}$ is the discrete Fourier transform matrix and $\mathcal{D}^\dagger$ denotes its Hermitian transpose, effectively implementing the inverse discrete Fourier transform operation. $\Psi_n$ is a diagonal matrix describing free space wave propagation by distance $z_n$ where the diagonal entries are raster-scanned from analytical values given by angular spectrum propagation kernel in the frequency domain, as defined in equation 1. $|−|^2$ is applied element-wise to compute intensity distribution from the propagated complex field.

To reconstruct the complex-valued elements of the object field from multi-plane intensity measurements, we pose the phase retrieval task as a nonlinear least-squares optimization problem. The objective is to minimize the discrepancy between the experimentally measured intensity images and the computed wavefield magnitudes derived from the forward model. In such a case, we define the cost function as,

$$ \phi(\mathbf{U}) = \frac{1}{2} \sum_{n=1}^{N} \left\| \mathbf{M}_n - |\mathcal{D}^\dagger \Psi_n \mathcal{D} \mathbf{U}|^2 \right\|_2^2 \quad (4) $$

where $\mathbf{M}_n \in \mathbb{R}^R$ is the vectorized form of $n^{th}$ measured intensity images and $||−||_2$ represents the Euclidean norm. Notably, the optimization problem in equation 4 is intrinsically non-convex, large-scale, and involves a modulus-squared nonlinearity, collectively making it computationally challenging to identify a direct
solution that effectively minimizes the corresponding cost function. To address this, we adopt an alternating minimization strategy, implemented through the iterative procedure. In such direction, we begin with an initial estimate of the complex field, $\mathbf{U}^0$, and then this estimate is subsequently updated at each iteration as follows,

$$ \mathbf{U}^{k+1} = \mathbf{U}^k + \alpha^k d^k \quad (5) $$

where $\alpha^k$ is the learning step size and $d^k$ represents the search direction. We employ a backtracking line search algorithm$^{50,56,63,64}$ to adaptively determine the step size, ensuring that the cost function decreases monotonically with each iteration. The next step involves evaluating a method for choosing the search direction, which is crucial for determining the convergence rate. A common approach is to use the first-order derivative (functional gradient) of the cost function. However, this implementation often suffers from slow convergence rates, a tendency to get trapped in local minima, and sensitivity to initialization, especially in ill-posed and non-convex phase retrieval problems$^{52}$. These issues can be alleviated by using a second-order optimization framework, which incorporates both the first-order derivative (gradient) and the second-order derivative (Hessian) of the cost function to improve convergence speed and reconstruction accuracy. Subsequently, similar to the quasi-Newton optimization approach$^{56,65}$, we evaluate the search direction by solving linear equations that incorporate the Hessian and the gradient of the cost function, as follows,

$$ \nabla^2_{(\mathbf{u}, \bar{\mathbf{u}})} \phi(\mathbf{U}, \bar{\mathbf{U}}) \left[ \frac{d}{d\bar{\mathbf{u}}} \right] + \nabla_{(\mathbf{u}, \bar{\mathbf{u}})} \phi(\mathbf{U}, \bar{\mathbf{U}}) = 0 \quad (6) $$

where $d$ is the search direction and $\bar{d}$ is the complex conjugate of $d$. Given the large size of the Hessian matrix in our problem, which scales with the square of the number of pixels, and computational complexity increases cubically with the number of pixels, directly solving equation 6 presents considerable computational challenges. In this context, Davidon-Fletcher-Powell, Broyden-Fletcher-Goldfarb-Shanno, and limited-memory Broyden-Fletcher-Goldfarb-Shanno methods have been introduced, which provide a crude approximation for the inverse of the Hessian matrix by leveraging estimates and gradients from previous iterations$^{56}$. However, these methods pose computational impracticalities, as they still necessitate the storage or computation of the Hessian matrix or its inverse, therefore, depend on the choice of memory size, i.e., the number of estimates and gradients retained in the iterations. To mitigate potential limitations arising from low memory, we employ an alternative approximation for the second-order derivative using the analytical expression for the gradient and Hessian of the cost function (see supplementary note S1), as follows,

$$ \nabla \phi(\mathbf{U}, \bar{\mathbf{U}}) = \sum_{n=1}^{N} \left[ \begin{matrix} -\mathcal{D}^{\dagger} \bar{\Psi}_{n} \mathcal{D} \mathrm{diag}(\mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} \mathbf{U}) \\ -\mathcal{D} \Psi_{n} \mathcal{D}^{\dagger} \mathrm{diag}(\mathcal{D} \bar{\Psi}_{n} \mathcal{D}^{\dagger} \bar{\mathbf{U}}) \end{matrix} \right] \left[ \mathbf{M}_{n} - |\mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} \mathbf{U}|^{2} \right] \\ \nabla^{2} \phi(\mathbf{U}, \bar{\mathbf{U}}) \approx \begin{bmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{bmatrix} \quad (7) $$

where,

$$ S_{11} = \mathcal{D}^{\dagger} \bar{\Psi}_{n} \mathcal{D} \mathrm{diag}(\mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} \mathbf{U}) \mathrm{diag}(\mathcal{D} \bar{\Psi}_{n} \mathcal{D}^{\dagger} \bar{\mathbf{U}}) \mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} $$

$$ S_{12} = \mathcal{D}^{\dagger} \bar{\Psi}_{n} \mathcal{D} \mathrm{diag}(\mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} \mathbf{U}) \mathrm{diag}(\mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} \mathbf{U}) \mathcal{D} \bar{\Psi}_{n} \mathcal{D}^{\dagger} $$

$$ S_{21} = \mathcal{D} \Psi_{n} \mathcal{D}^{\dagger} \mathrm{diag}(\mathcal{D} \bar{\Psi}_{n} \mathcal{D}^{\dagger} \bar{\mathbf{U}}) \mathrm{diag}(\mathcal{D} \bar{\Psi}_{n} \mathcal{D}^{\dagger} \bar{\mathbf{U}}) \mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} $$

$$ S_{22} = \mathcal{D} \Psi_{n} \mathcal{D}^{\dagger} \mathrm{diag}(\mathcal{D} \bar{\Psi}_{n} \mathcal{D}^{\dagger} \bar{\mathbf{U}}) \mathrm{diag}(\mathcal{D}^{\dagger} \Psi_{n} \mathcal{D} \mathbf{U}) \mathcal{D} \bar{\Psi}_{n} \mathcal{D}^{\dagger} $$

Equation 7 provides a direct analytical expression for the gradient and Hessian of the cost function, facilitating the solution of the linear equations (equation 6) via an iterative conjugate gradient method$^{53}$. It is worth noting that computing solutions for $d^k$ using an analytical expression incurs no additional memory overhead. Furthermore, we introduce a subroutine to evaluate $d^k$ based on the current estimate $\mathbf{U}^k$, which is iteratively updated in the object field reconstruction step (equation 5). The accuracy of this subroutine, however, depends inherently on the number of iterations within the inner solver, thereby introducing a trade-off between computational efficiency and convergence rate. Despite this, our search direction evaluation benefits from fast matrix-vector multiplications implemented using the Fast Fourier Transform,
achieving a computational complexity of $O(N\log N)$. Such enhancement significantly accelerates
convergence and reduces overall runtime.

## Results and Discussion

### Numerical simulations and performance analysis

We evaluate the efficiency of our second-order phase retrieval framework by conducting a series of numerical simulations, as illustrated in Figure 2. In this context, we simulated an imaging target comprising both amplitude (Figure 2a) and phase (Figures 2b and 2c) components. In our simulations, designed to closely mimic realistic experimental conditions, we employ a source wavelength of 85.7 µm (equivalent to 3.5 THz frequency), the imaging field-of-view is discretized into 100 × 100 pixels, with a transverse spatial sampling interval of 50 µm, and spanning the sample illuminated area of 5 × 5 mm². In the lens-less imaging settings, the THz camera is placed at 2.5 mm from the object plane, and a sequence of 11 diffraction patterns is numerically simulated at equidistant axial steps of 0.3 mm. Conversely, for defocusing-based imaging configuration, we numerically simulated 11 intensity measurements at incremental defocus planes with an axial spacing of 0.2 mm symmetrically around the focal plane. Furthermore, to assess the robustness of our numerical framework under realistic experimental conditions, we incorporated additive white-Gaussian noise into the simulations while emulating all the intensity patterns, corresponding to a SNR of 20 dB, consistent with performance metrics typically observed in state-of-the-art THz detection microbolometric cameras. Interestingly, the amplitude and phase profiles are numerically reconstructed in both a lens-less imaging scheme (Figures 2d-f) and a defocusing-based imaging configuration (Figures 2g-i), demonstrating

<figure><img src="image_3.png" /><figcaption>**Figure 2: Performance analysis of a second-order optimization-based phase retrieval scheme. (a)**</figcaption></figure>

Amplitude, (b) phase, and (c) slowly varying phase of the simulated true imaging sample. (d-f) Reconstructed amplitude and phase distributions in a lens-less imaging configuration. (g-i) Reconstructed amplitude and phase in a defocusing imaging configuration. (j) Vertical line cut comparing spatial phase variation in the test image (dotted black), lens-less (red), and defocusing imaging configurations (blue). (k) Variation of the phase structural similarity index with the number of intensity measurements required for lens-less (red) and defocusing (blue) imaging configurations. (l) Variation of the phase structural similarity index as a function of SNR detection error for 11 emulated intensity patterns in lens-less (red) and defocused (blue) imaging configurations. The 5 × 5 mm² sample illumination area is spatially sampled at 50 µm resolution, corresponding to 100 × 100 pixels.
the versatility of the proposed phase retrieval method. Additionally, to rigorously assess the robustness of our technique, we extended our analysis to retrieve complex spatial phase profiles characterized by slowly varying spatial variations. In our simulations, this scenario is emulated by modeling a spatial phase distribution as a circular Gaussian random variable with a $\mu = 0, \sigma = \frac{1}{2}$. To introduce controllable spatial features, we applied a Gaussian filter with a standard deviation of 200µm, effectively mimicking a physically realizable passive phase object (Figure 2c). The corresponding phase reconstructions for such randomized phase target are shown in Figures 2f and 2i. Our approach successfully recovers the spatial phase morphology in both configurations. To quantitatively validate reconstruction fidelity, we performed line profile comparisons along the vertical (y-axis) direction, as marked by the dashed lines in Figures 2c, 2f, and 2i. These comparisons, summarized in Figure 2j, show excellent agreement between the reconstructed and true phase profiles and confirm the reliability and precision of our second-order retrieval method in complex phase imaging scenarios.

Furthermore, to evaluate the performance of our proposed method, we analyzed the minimum number of intensity patterns required for accurate phase reconstruction under both imaging configurations. For a quantitative assessment of reconstruction quality, we employed the Structural Similarity Index (SSIM)⁶⁶ metric to compare the reconstructed phase profiles with the corresponding test image. As shown in Figure 2k, the SSIM increases with the number of intensity patterns, indicating improved reconstruction fidelity. Notably, our approach achieves maximal structural similarity with as few as 11 intensity patterns, highlighting its efficiency in terms of data acquisition. In addition, we investigated the robustness of the proposed method under varying detection noise levels. Specifically, we introduced additive white Gaussian noise into the 11 simulated intensity measurements to emulate realistic experimental conditions and evaluate the robustness of the reconstruction across different detection SNR levels. Figure 2l illustrates the SSIM variation as a function of SNR detection error. The results demonstrate that the algorithm maintains stable performance even under low SNR conditions, achieving reliable reconstruction at an SNR threshold as low as 15 dB. Such robustness demonstrates the practical applicability and resilience of our phase retrieval framework in noisy measurement conditions.

## Experimental verification

In the experimental configurations shown in Figures 1a and 1b, we employ a CW QCL source from Lytid, operating at a monochromatic frequency of 3.5 THz (equivalent to a wavelength of 85.7 µm), which emits a divergent beam with an output power of ~13 mW (although we use 70% of source power in the experiment). To achieve beam collimation, we utilize a beam-collimating optical system consisting of three off-axis parabolic mirrors with focal lengths of 15.24 cm, 22.86 cm, and 7.62 cm, in conjunction with a pinhole. Conversely, the test sample is illuminated with a spatially coherent beam of diameter 0.68 cm, ensuring that the transverse coherence area is larger than the size of the features examined in the sample. The transmitted field from the sample is recorded across multiple axial planes using an uncooled microbolometric detector array (INO MICROXCAM-384i), which provides 384 × 288 pixels images, with a pixel pitch of 35 µm, and a noise-equivalent power (NEP) of ~0.15 pW/√Hz⁶⁷.

The first experimental configuration employs a lens-less scheme (Figure 1a) where the microbolometer, initially positioned 3.5 mm from the object, is incrementally translated along the optical axis in 0.2 mm steps to record a series of 11 diffraction patterns corresponding to different propagation distances. Each recorded intensity frame is trimmed to a 250 × 250 pixels region of interest, encompassing a field-of-view of 8.75 × 8.75 mm² with a pixel pitch of 35 µm. In the transmission mode of the lens-less imaging configuration, the wings of a *European hornet* and *Archaeophya adamsi dragonfly* (Figure 3a and 3f) are used as imaging samples due to their complex topographical features, which fall within the lateral resolution range and provide an ideal testbed for evaluating the performance of phase reconstruction. Figure 3b and 3g, represent a set of 11 diffraction patterns at collected at multiple planes for the test objects shown in Figure 3a and 3f, respectively. The spatial phase profiles reconstructed using our proposed second-order optimization-based algorithm are shown in Figures 3c and 3h, revealing detailed structural features with high fidelity. For comparison, Figures 3d and 3i represent the corresponding phase reconstructions obtained
<figure><img src="image_4.png" /><figcaption>**Figure 3:** THz lens-less phase imaging using multi-plane intensity measurements. (a, f) Optical images of the *European hornet* and *Archaeophya adamsi* insect wings as imaging samples. (b, g) recorded multi-plane diffraction patterns for the test object in (a) and (f), respectively. (c, h) Reconstructed spatial phase distribution obtained using the proposed second-order optimization approach. (d, i) Phase reconstructions obtained from a conventional multi-plane phase retrieval method. (e, j) The variation of the cost function (mean square error between the measured and predicted intensity) with respect to the number of iterations for the test sample presented in (a) and (f), respectively. The 8.75 × 8.75 mm² sample illumination area is spatially sampled at 35 µm resolution, corresponding to 250 × 250 pixels.</figcaption></figure>

via a standard multiple-plane phase retrieval method described in the reference⁴³,⁶⁸. The contrast between our approach and multiple-plane phase retrieval is evident in the phase reconstruction results. As shown in Figures 3d and 3i, the standard multiple plane phase retrieval method shows phase inaccuracies, particularly in the reconstruction of low spatial frequency components. The quantitative assessment of convergence, illustrated through the evolution of the cost function over iterations, is shown in Figures 3e and 3j. The convergence metric (cost function) illustrates that the second-order optimization approach demonstrates a more rapid convergence compared to the multi-plane phase retrieval technique. Specifically, with our
second-order optimization, the cost function decreases significantly faster and ultimately attains a notably lower steady-state value, highlighting better performance in minimizing the reconstruction error. Additionally, our proposed approach achieves convergence with significantly fewer iterations, emphasizing its computational efficiency and reduced runtime relative to the iterative multi-plane phase retrieval strategy.

In the second configuration (Figure 1b), a 4f imaging system is implemented into the detection path using a pair of identical silicon plano-convex lenses ($f = 3.5$ cm), with the detector positioned at the focal plane.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** THz phase retrieval from a stack of defocusing intensity measurements. (a, f) Optical images of imaging samples. (b, g) Measured stack of defocused intensity patterns for the test object shown in (a) and (f), respectively. (c, h) Reconstructed spatial phase distribution obtained using the proposed second-order optimization approach. (d, i) Phase reconstructions obtained from conventional multi-plane phase retrieval method. (e, j) The variation of the cost function (mean square error between the measured and predicted intensity) with respect to the number of iterations for the test sample presented in (a) and (f), respectively. The $7 \times 7 \text{ mm}^2$ sample illumination area is spatially sampled at 35 µm resolution, corresponding to $200 \times 200$ pixels.</figcaption></figure>
In such a case, the defocused wavefront intensities are acquired at multiple planes symmetrically distributed around the focal plane, with a constant axial step size of 0.2 mm. Such deliberate defocusing introduces spatial interference fringes at the image plane and compromises image sharpness due to defocus-induced blurring. Further, the intensity images recorded at each plane are trimmed to 200 × 200 pixels, corresponding to a 7 × 7 mm² region of interest, with a spatial sampling of 35 µm. In our experiment, we used a *polypropylene plastic 'Arrow' and a 'gear-shaped lock washer'* as test samples, as shown in Figures 4a and 4f, respectively. Figures 4b and 4g show recorded defocused intensity patterns for the test objects, as shown in Figures 4a and 4f. The reconstructed spatial phase distributions obtained through our proposed approach are shown in Figures 4c and 4h. For comparison, we also retrieved the spatial phase profiles using the iterative multi-plane phase retrieval algorithm<sup>69–71</sup>, shown in Figures 4d and 4i. Notably, the phase reconstructions obtained via our technique exhibit better fidelity than the iterative multi-plane phase retrieval technique. The fidelity of such reconstruction can be quantitatively assessed by examining the convergence behavior of the cost function, as shown in Figures 4e and 4j. Our proposed second-order optimization technique demonstrates significantly faster convergence compared to the multi-plane phase retrieval method. Notably, our approach demonstrates a steep decline in the cost function, reaching a lower constant value than that of the multi-plane phase retrieval method, which indicates improved accuracy in minimizing the cost function.

Both experimental configurations are conceptually aligned in acquiring multiple intensity measurements across varying propagation distances; however, they diverge notably in their practical implementation. The lens-less arrangement offers a minimalistic architecture, requiring fewer optical components and enabling a more compact and straightforward setup. This configuration facilitates direct acquisition of diffraction patterns and simplifies the data acquisition workflow. In contrast, the 4f telecentric system is well-suited for imaging large specimens within the THz spectral range, particularly when limited field-of-view detectors are utilized. In this approach, the axial translation steps of the detector can be systematically determined by geometric scaling transformation as $z'_{\pm} = z_{\pm}(f_1/f_2)$, where $f_1$ and $f_2$ are focal lengths of lenses in the 4f telecentric imaging system. Despite its utility, the 4f configuration imposes higher energy demands due to the necessity for a spatially extended THz beam and the inclusion of a large-aperture lens. While fabricating such lenses with high surface quality and low roughness poses significant challenges in the visible spectrum, these constraints are comparatively relaxed in the THz domain. Ultimately, employing both configurations in our study demonstrates the adaptability and broad applicability of the proposed computational reconstruction framework under varying optical and imaging system-level constraints. Additionally, regarding experiment timeframe, the MATLAB-integrated automatic data acquisition module facilitates recording of 11 intensity patterns for both configurations in ~45 seconds. The subsequent second-order optimization procedure, involving 50 iterations and analytical computation of the gradient and Hessian via linear equation solving, is completed in ~10 seconds on a MacBook Pro (M4, GPU-accelerated 48 GB RAM).

To demonstrate the practical applicability of our THz phase imaging approach, we investigate its capacity to spatially resolve chemical heterogeneity, specifically regarding water content distribution within leaf samples. In this context, we employ our phase reconstruction technique to image a leaf (as shown in Figure 5a) and quantify the internal water distribution by analyzing the spatial phase variations. While THz absorption due to water strongly affects amplitude and phase, phase imaging is often more sensitive to subtle structural and compositional changes since it directly measures refractive index variations and optical path differences. On the other hand, direct single-shot intensity imaging provides semi-quantitative information because amplitude reduction can often be masked by scattering, reflections, and surface roughness, resulting in degraded image contrast. Therefore, the reconstructed phase is a reliable proxy for mapping its presence<sup>72</sup>. To this end, we compare the phase morphology of a semi-fresh leaf with that of the same specimen after undergoing 72 hours of dehydration. In the leaf, the region of interest spans 8.75 × 8.75 mm², with phase maps reconstructed at a resolution of 250 × 250 pixels, corresponding to a spatial sampling of 35 µm. A set of defocused intensity patterns for fresh and dried leaf is shown in Figures 5b and 5c, respectively. Figures 5d and 5e illustrate the retrieved phase morphology of the fresh and desiccated
<figure><img src="image_6.png" /><figcaption>**Figure 5:** THz phase imaging of a leaf. (a) Optical image of a leaf. (b, c) recorded set of defocused intensity patterns for fresh and dried leaf, respectively. (d) Reconstructed spatial phase map of fresh leaf. (e) Reconstructed spatial phase distribution of dried leaf. The 8.75 × 8.75 mm² sample illumination area is spatially sampled at 35 µm resolution, corresponding to 250 × 250 pixels.</figcaption></figure>

leaves, respectively. It is worth noting that despite the macroscopic appearance of the leaf remaining largely unaltered, the phase reconstructions show substantial internal structural changes. Remarkably, the observed differences demonstrate a significant redistribution of water content, particularly evident by increased contrast along the vascular bundles, indicating that much of the water has evaporated from the surrounding tissue during the monitored time frame.

## Conclusion

In this work, we have presented and experimentally validated a non-interferometric framework for phase recovery of samples in the THz spectral range, leveraging a robust nonlinear optimization scheme. Our methodology introduces two practical imaging configurations, a lens-less and a 4f imaging configuration,
that enable the acquisition of multi-plane intensity patterns without requiring significant hardware modifications. We demonstrated the practical applicability of our phase retrieval technique for synthesizing the internal structural composition of the sample. As a relevant case example, we reconstructed spatial phase profiles of both natural and dehydrated leaves, allowing for an estimation of hydration levels: a promising application in biosensing. It is worth noting that our imaging approach fundamentally differs from TDS-based imaging, where both field amplitude and phase are directly accessible without the need for computational phase retrieval. Whereas, our technique, based on CW illumination and intensity-only detection, is specifically designed for scenarios where direct phase access is not feasible. At the core of our approach is a second-order modified Gauss-Newton optimization algorithm, utilizing analytically derived expressions for both the gradient and Hessian of the cost function. This formulation effectively overcomes the limitations of previous first-order methods, which often suffer from slow convergence rates and excessive memory consumption. Comparative evaluations demonstrate that our method surpasses classical algorithms, such as standard multi-plane phase retrieval, both in terms of convergence speed and reconstruction quality. Notably, our approach exhibits remarkable robustness against noise and maintains reliable performance even under low-detection SNR conditions. As a relevant case example in our numerical analysis, spatial phase reconstruction achieves an SSIM of ~1 even at a detection SNR of 15 dB. Additionally, our methodology demonstrated computational efficiency: for instance, data acquisition for multiple planes can be completed in approximately 45 seconds, while 50 iterations to approximate the Hessian of the cost function by solving the Hessian linear equation take ~10 seconds, resulting in a total phase reconstruction time of ~1 minute. While our current implementation acquires intensity patterns at uniformly spaced axial positions to ensure redundancy and diversity in the data, the flexibility in plane selection remains an advantage. Fine spatial phase details (high spatial frequency components) are optimally captured at smaller plane separations, whereas broader phase variations (low spatial frequency components) are most effectively retrieved at larger separations, offering a tunable framework adaptable to diverse imaging requirements. The initial guess is pivotal in ensuring efficient convergence in second-order Gauss-Newton-type optimization schemes. A promising direction to significantly reduce the number of iterations required and accelerate the convergence of the proposed phase retrieval framework involves leveraging the initial estimate obtained from the multiplane phase retrieval solution. As a future direction, an alternative algorithmic pathway worth exploring is the stochastic gradient descent, which could offer improved memory efficiency, making it particularly advantageous for large-scale or high-throughput imaging scenarios. Our THz phase imaging paradigm not only provides innovative avenues for non-destructive investigation of internal sample structures but also opens promising opportunities across various fields, including biomedical imaging, chemical analysis, and materials science.

**Supplementary Material**

Supplementary note 1: Derivation of first-order and second-order derivatives of the cost function with CR calculus

**Funding**

Agence Nationale de la Recherche (ANR-22CE42-0005-HYPSTER, ANR 22-PEEL-0003-Comptera).

**Disclosures**

The authors declare no conflicts of interest

**Data Availability**

The datasets for all figures are freely accessible at [figshare]().

**Author contributions**

V. K. conceived the idea. All authors were engaged in the discussion regarding the basic concept of the manuscript and its experimental implementation for THz waves. Experimental investigators: V.K. (lead), P.Mu. (equal), and F.F. (supporting). Formal analysis and data curation: V. K. (lead) and P. Mu. (equal).
Investigation and results interpretation: V.K. (lead), P.Mu. (equal), K.K. (equal), S.G. (equal), and P.Mo. (equal). Manuscript drafting: V.K. (lead), P.Mu. (supporting), K.K. (supporting), S.G. (supporting), and P.Mo. (supporting). Project supervision: P. Mo. (equal) and S. G. (equal).

References

¹ D.M. Mittleman, “Twenty years of terahertz imaging [Invited],” Opt. Express **26**(8), 9417 (2018).

² E. Castro-Camus, M. Koch, and D.M. Mittleman, “Recent advances in terahertz imaging: 1999 to 2021,” Appl. Phys. B **128**(1), 12 (2021).

³ G. Valušis, A. Lisauskas, H. Yuan, W. Knap, and H.G. Roskos, “Roadmap of Terahertz Imaging 2021,” Sensors **21**(12), 4092 (2021).

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

⁵ V. Anitha, A. Beohar, and A. Nella, “THz Imaging Technology Trends and Wide Variety of Applications: a Detailed Survey,” Plasmonics **18**(2), 441–483 (2023).

⁶ K. Khare, M. Butola, and S. Rajora, *Fourier Optics and Computational Imaging*, 2nd ed. (Springer-Nature, 2023).

⁷ M. Wan, J.J. Healy, and J.T. Sheridan, “Terahertz phase imaging and biomedical applications,” Opt. Laser Technol. **122**, 105859 (2020).

⁸ J. Neu, and C.A. Schmuttenmaer, “Tutorial: An introduction to terahertz time domain spectroscopy (THz-TDS),” J. Appl. Phys. **124**(23), 231101 (2018).

⁹ V. Jelic, S. Adams, M. Hassan, K. Cleland-Host, S.E. Ammerman, and T.L. Cocker, “Atomic-scale terahertz time-domain spectroscopy,” Nat. Photonics **18**(9), 898–904 (2024).

¹⁰ Y. Cao, P. Huang, P. Huang, J. Chen, W. Ge, D. Hou, G. Zhang, and G. Zhang, “Qualitative and quantitative detection of liver injury with terahertz time-domain spectroscopy,” Biomed. Opt. Express **11**(2), 982–993 (2020).

¹¹ A. Chen, A. Chen, O.B. Osman, O.B. Osman, Z.B. Harris, A. Abazri, R. Honkanen, and M.H. Arbab, “Investigation of water diffusion dynamics in corneal phantoms using terahertz time-domain spectroscopy,” Biomed. Opt. Express **11**(3), 1284–1297 (2020).

¹² Q. Cassar, S. Caravera, G. MacGrogan, T. Bücher, P. Hillger, U. Pfeiffer, T. Zimmer, J.-P. Guillet, and P. Mounaix, “Terahertz refractive index-based morphological dilation for breast carcinoma delineation,” Sci. Rep. **11**(1), 6457 (2021).

¹³ 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**(3), 687–694 (2020).

¹⁴ L.L. Hale, T. Siday, and O. Mitrofanov, “Near-field imaging and spectroscopy of terahertz resonators and metasurfaces [Invited],” Opt. Mater. Express **13**(11), 3068–3086 (2023).

¹⁵ V. Kumar, P. Mukherjee, L. Valzania, A. Badon, P. Mounaix, and S. Gigan, “Fourier Synthetic Aperture-based Time-resolved Terahertz Imaging,” Photonics Res. **13**(2), 407–416 (2024).

¹⁶ 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 Sci. Appl. **8**(1), 55 (2019).

¹⁷ 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**(5), 521 (2020).
18 L. Olivieri, J.S.T. Gongora, L. Peters, V. Cecconi, A. Cutrona, J. Tunesi, R. Tucker, A. Pasquazi, and M. Peccianti, “Hyperspectral terahertz microscopy via nonlinear ghost imaging,” *Optica* **7**(2), 186–191 (2020).

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

20 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**(8), 2634–2642 (2022).

21 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* **11**(2), 362–368 (2024).

22 V. Kumar, V. Cecconi, A. Cutrona, L. Peters, L. Olivieri, J.S. Totero Gongora, A. Pasquazi, and M. Peccianti, “Terahertz microscopy through complex media,” *Sci. Rep.* **15**(1), 11706 (2025).

23 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.* **36**(11), 1993–1995 (2011).

24 P. Zolliker, and E. Hack, “THz holography in reflection using a high resolution microbolometer array,” *Opt. Express* **23**(9), 10957–10967 (2015).

25 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**(1), 13566 (2015).

26 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* **26**(20), 25805–25813 (2018).

27 K. Xue, Q. Li, Y.-D. Li, and Q. Wang, “Continuous-wave terahertz in-line digital holography,” *Opt. Lett.* **37**(15), 3228–3230 (2012).

28 L. Rong, T. Latychevskaja, 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**(1), 8445 (2015).

29 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,” *Opt. Commun.* **426**, 612–622 (2018).

30 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* **27**(2), 702–713 (2019).

31 J.T. Sheridan, R.K. Kostuk, A.F. Gil, Y. Wang, W. Lu, H. Zhong, Y. Tomita, C. Neipp, J. Francés, S. Gallego, I. Pascual, V. Marinova, S.-H. Lin, K.-Y. Hsu, F. Bruder, S. Hansen, C. Manecke, R. Meisenheimer, C. Rewitz, T. Rölle, S. Odinokov, O. Matoba, M. Kumar, X. Quan, Y. Awatsuji, P.W. Wachulak, A.V. Gorelaya, A.A. Sevryugin, E.V. Shalymov, V. Yu Venediktov, R. Chmelik, M.A. Ferrara, G. Coppola, A. Márquez, A. Beléndez, W. Yang, R. Yuste, A. Bianco, A. Zanutta, C. Falldorf, J.J. Healy, X. Fan, B.M. Hennelly, I. Zhurminsky, M. Schnieper, R. Ferrini, S. Fricke, G. Situ, H. Wang, A.S. Abdurashitov, V.V. Tuchin, N.V. Petrov, T. Nomura, D.R. Morim, and K. Saravanamuttu, “Roadmap on holography,” J. Opt. **22**(12), 123002 (2020).
32 R. Ivaškevičiūtė-Povilauskienė, I. Grigelionis, A. Siemion, D. Jokubauskis, K. Ikamas, A. Lisauskas, L. Minkevičius, and G. Valušis, “Multifrequency Digital Terahertz Holography Within 1.39–4.25 THz Range,” IEEE Trans. Terahertz Sci. Technol. 14(5), 568–578 (2024).

33 M.S. Heimbeck, M.K. Kim, D.A. Gregory, and H.O. Everitt, “Terahertz digital holography using angular spectrum and dual wavelength reconstruction methods,” Opt. Express 19(10), 9192–9200 (2011).

34 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 24(18), 21134–21146 (2016).

35 L. Valzania, P. Zolliker, and E. Hack, “Coherent reconstruction of a textile and a hidden object with terahertz radiation,” Optica 6(4), 518–523 (2019).

36 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. 55(3), A43–A48 (2016).

37 S. Krimi, J. Klier, J. Jonuscheit, G. von Freymann, R. Urbansky, and R. Beigang, “Highly accurate thickness measurement of multi-layered automotive paints using terahertz technology,” Appl. Phys. Lett. 109(2), 021105 (2016).

38 M. Xiang, H. Yuan, L. Wang, K. Zhou, and H.G. Roskos, “Amplitude/Phase Retrieval for Terahertz Holography With Supervised and Unsupervised Physics-Informed Deep Learning,” IEEE Trans. Terahertz Sci. Technol. 14(2), 208–215 (2024).

39 L. Valzania, T. Feurer, P. Zolliker, and E. Hack, “Terahertz ptychography,” Opt. Lett. 43(3), 543–546 (2018).

40 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 27(2), 938–950 (2019).

41 V. Kumar, P. Mukherjee, F. Fauquet, A. Badon, P. Mounaix, and S. Gigan, “Terahertz Fourier Ptychography for Complex Media Imaging,” in 2024 49th Int. Conf. Infrared Millim. Terahertz Waves IRMMW-THz, (2024), pp. 1–3.

42 L. Peng, C. Luo, B. Deng, H. Wang, S. Chen, and J. Dong, “Phaseless Terahertz Coded-Aperture Imaging for Sparse Target Based on Phase Retrieval Algorithm,” Sensors 19(21), 4617 (2019).

43 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. 45(15), 4168–4171 (2020).

44 C. Xing, F. Qi, and S. Guo, “Enhanced Terahertz Phase Retrieval Imaging by Unequal Spaced Measurement,” Sensors 22(10), 3816 (2022).

45 A. Chopard, E. Tsiplakova, N. Balbekin, O. Smolyanskaya, J.-B. Perraud, J.-P. Guillet, N.V. Petrov, and P. Mounaix, “Single-scan multiplane phase retrieval with a radiation of terahertz quantum cascade laser,” Appl. Phys. B 128(3), 63 (2022).

46 J.R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982).

47 M.R. Teague, “Deterministic phase retrieval: a Green’s function solution,” JOSA 73(11), 1434–1441 (1983).

48 N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. 49(1), 6–10 (1984).

49 G.R. Brady, and J.R. Fienup, “Nonlinear optimization algorithm for retrieving the full complex pupil function,” Opt. Express 14(2), 474–486 (2006).
50 S. Boyd, and L. Vandenberghe, *Convex Optimization* (Cambridge University Press, 2004).

51 D. Kim, and J.A. Fessler, “Optimized first-order methods for smooth convex minimization,” Math. Program. **159**(1), 81–107 (2016).

52 J. Dong, L. Valzania, A. Maillard, T. Pham, S. Gigan, and M. Unser, “Phase Retrieval: From Computational Imaging to Machine Learning: A tutorial,” IEEE Signal Process. Mag. **40**(1), 45–57 (2023).

53 J.L. Nazareth, “Conjugate gradient method,” WIREs Comput. Stat. **1**(3), 348–353 (2009).

54 G. Qu, and N. Li, “Accelerated distributed Nesterov Gradient Descent for convex and smooth functions,” in *2017 IEEE 56th Annu. Conf. Decis. Control CDC*, (2017), pp. 2260–2267.

55 G. Qu, and N. Li, “Accelerated Distributed Nesterov Gradient Descent,” IEEE Trans. Autom. Control **65**(6), 2566–2581 (2020).

56 J. Nocedal, and S. Wright J., *Numerical Optimization* (Springer New York, 2006).

57 D.C. Liu, and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Program. **45**(1), 503–528 (1989).

58 M. Guizar-Sicairos, and J.R. Fienup, “Phase retrieval with transverse translation diversity: a nonlinear optimization approach,” Opt. Express **16**(10), 7264–7278 (2008).

59 J. Zhong, L. Tian, P. Varma, and L. Waller, “Nonlinear Optimization Algorithm for Partially Coherent Phase Retrieval and Source Recovery,” IEEE Trans. Comput. Imaging **2**(3), 310–322 (2016).

60 T. Wu, J. Dong, X. Shao, and S. Gigan, “Imaging through a thin scattering layer and jointly retrieving the point-spread-function using phase-diversity,” Opt. Express **25**(22), 27182–27194 (2017).

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

62 K. Matsushima, and T. Shimobaba, “Band-Limited Angular Spectrum Method for Numerical Simulation of Free-Space Propagation in Far and Near Fields,” Opt. Express **17**(22), 19662–19673 (2009).

63 L. Armijo, “Minimization of functions having Lipschitz continuous first partial derivatives,” Pac. J. Math. **16**(1), 1–3 (1966).

64 D.P. Bertsekas, and J.N. Tsitsiklis, “Gradient Convergence in Gradient methods with Errors,” SIAM J. Optim. **10**(3), 627–642 (2000).

65 F.J. Aragón, M.A. Goberna, M.A. López, and M.M.L. Rodríguez, *Nonlinear Optimization* (Springer International Publishing, Cham, 2019).

66 Z. Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. **13**(4), 600–612 (2004).

67 M. Terroux, P. Talbot, F. Généreux, L. Marchese, E.-H. Oulachgar, and A. Bergeron, “NEP characterization and analysis method for THz imaging devices,” in *Passive Act. Millim.-Wave Imaging XXIV*, (SPIE, 2021), pp. 146–151.

68 E.G. Tsiplakova, J.-B. Perraud, J.-P. Guillet, P. Mounaix, and N.V. Petrov, “Terahertz diffractive imaging with saturated data inpainting: publisher’s note,” Opt. Lett. **50**(6), 1844–1844 (2025).

69 R. Gerchberg, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik **35**, 237 (1972).

70 L.J. Allen, and M.P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. **199**(1), 65–75 (2001).
<footnote n="71">⁷¹ P. Memmolo, L. Miccio, F. Merola, A. Paciello, V. Embrione, S. Fusco, P. Ferraro, and P. Antonio Netti, “Investigation on specific solutions of Gerchberg–Saxton algorithm,” Opt. Lasers Eng. **52**, 206–211 (2014).</footnote>

<footnote n="72">⁷² B.B. Hu, and M.C. Nuss, “Imaging with terahertz waves,” Opt. Lett. **20**(16), 1716–1718 (1995).</footnote>
# Supplementary materials

## Computational Terahertz Phase Imaging using a Sequence of Multi-Plane Intensity Measurements

Vivek Kumar<fnref n="1" />,<fnref n="2" />, Pitambar Mukherjee<fnref n="2" />, Frédéric Fauquet<fnref n="2" />, Kedar Khare<fnref n="3" />, Sylvain Gigan<fnref n="1" />, Patrick Mounaix<fnref n="2" />

## Supplementary note S1

We reformulate the general representation of the cost function for the phase recovery problem as defined in equation 4:

$$ \phi(\mathbf{U}) = \frac{1}{2} \sum_{n=1}^{N} (\mathbf{M}_n - |T_n \mathbf{U}|^2)^{\dagger} (\mathbf{M}_n - |T_n \mathbf{U}|^2) = \sum_{n=1}^{N} \mathcal{L}_n(\mathbf{U}, \bar{\mathbf{U}}) \quad (1) $$

Here, $T_n = \mathcal{D}^\dagger \Psi_n \mathcal{D}$ denotes the transfer function of the free space, $\mathcal{L}_n(\mathbf{U}, \bar{\mathbf{U}})$ term represents error for $n^{th}$ intensity pattern, and it can be defined as,

$$ \mathcal{L}_n(\mathbf{U}, \bar{\mathbf{U}}) = \frac{1}{2} \ell_n(\mathbf{U}, \bar{\mathbf{U}})^{\dagger} \ell_n(\mathbf{U}, \bar{\mathbf{U}}) = \frac{1}{2} \sum_{r=1}^{R} \ell_{nr}(\mathbf{U}, \bar{\mathbf{U}}) \bar{\ell}_{nr}(\mathbf{U}, \bar{\mathbf{U}}) \quad (2) $$

where $\ell_n(\mathbf{U}, \bar{\mathbf{U}}) = \mathbf{M}_n - |T_n \mathbf{U}|^2$ is a vector representing the difference between measurement and prediction for $n^{th}$ intensity pattern, $R$ is the length of $\ell_n(\mathbf{U}, \bar{\mathbf{U}})$ vector, and $\ell_{nr}(\mathbf{U}, \bar{\mathbf{U}})$ is $r^{th}$ entry in the vector. Further, the first and second-order derivatives for each component $\mathcal{L}_n(\mathbf{U}, \bar{\mathbf{U}})$ can be derived using the chain rule and written as,

$$ \nabla \mathcal{L}_n(\mathbf{U}, \bar{\mathbf{U}}) = \mathbf{J}_n(\mathbf{U}, \bar{\mathbf{U}})^{\dagger} \ell_n(\mathbf{U}, \bar{\mathbf{U}}) \\ \nabla^2 \mathcal{L}_n(\mathbf{U}, \bar{\mathbf{U}}) = \mathbf{J}_n(\mathbf{U}, \bar{\mathbf{U}})^{\dagger} \mathbf{J}_n(\mathbf{U}, \bar{\mathbf{U}}) + \sum_{r=1}^{R} \ell_{nr}(\mathbf{U}, \bar{\mathbf{U}}) \ell''_{nr}(\mathbf{U}, \bar{\mathbf{U}}) \quad (3) $$

$$ \nabla^2 \mathcal{L}_n(\mathbf{U}, \bar{\mathbf{U}}) = \mathbf{J}_n(\mathbf{U}, \bar{\mathbf{U}})^{\dagger} \mathbf{J}_n(\mathbf{U}, \bar{\mathbf{U}}) + \sum_{r=1}^{R} \ell_{nr}(\mathbf{U}, \bar{\mathbf{U}}) \ell''_{nr}(\mathbf{U}, \bar{\mathbf{U}}) \quad (3) $$

where $\mathbf{J}_n(\mathbf{U}, \bar{\mathbf{U}})^{\dagger}$ is the Jacobian of the function vector $\ell_n(\mathbf{U}, \bar{\mathbf{U}})$ and can be obtained from $\mathbb{C}\mathbb{R}$ calculus theory [1] as,

$$ \mathbf{J}_n(\mathbf{U}, \bar{\mathbf{U}}) = [-\mathrm{diag}(\bar{T}_n \bar{\mathbf{U}})T_n - \mathrm{diag}(T_n \mathbf{U})\bar{T}_n] \quad (4) $$

We use the Gauss-Newton method [2] to approximate the Hessian matrix using only the Jacobian matrix, while omitting the second term in equation 9. Further, by substituting the value for $T_n$, the final analytical expressions for the gradient and Hessian of the cost function can be written as,

$$ \nabla \phi(\mathbf{U}, \bar{\mathbf{U}}) = \sum_{n=1}^{N} \left[ \begin{matrix} -\mathcal{D}^{\dagger} \bar{\Psi}_n \mathcal{D} \mathrm{diag}(\mathcal{D}^{\dagger} \Psi_n \mathcal{D} \mathbf{U}) \\ -\mathcal{D} \Psi_n \mathcal{D}^{\dagger} \mathrm{diag}(\mathcal{D} \bar{\Psi}_n \mathcal{D}^{\dagger} \bar{\mathbf{U}}) \end{matrix} \right] \left[ \mathbf{M}_n - |\mathcal{D}^{\dagger} \Psi_n \mathcal{D} \mathbf{U}|^2 \right] \\ \nabla^2 \phi(\mathbf{U}, \bar{\mathbf{U}}) \approx \begin{bmatrix} S_{11} & S_{12} \\ S_{21} & S_{22} \end{bmatrix} \quad (5) $$

where,

<footnote n="1">¹Laboratoire Kastler Brossel, ENS-Université PSL, CNRS, Sorbonne Université, Collège 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">³Optics and Photonics Centre, Indian Institute of Technology Delhi, New Delhi 110016, India.</footnote>
$$S_{11} = \mathcal{D}^{\dagger}\bar{\Psi}_{n}\mathcal{D} \operatorname{diag}(\mathcal{D}^{\dagger}\Psi_{n}\mathcal{D}\mathbf{U})\operatorname{diag}(\mathcal{D}\bar{\Psi}_{n}\mathcal{D}^{\dagger}\bar{\mathbf{U}})\mathcal{D}^{\dagger}\Psi_{n}\mathcal{D}$$

$$S_{12} = \mathcal{D}^{\dagger}\bar{\Psi}_{n}\mathcal{D} \operatorname{diag}(\mathcal{D}^{\dagger}\Psi_{n}\mathcal{D}\mathbf{U})\operatorname{diag}(\mathcal{D}^{\dagger}\Psi_{n}\mathcal{D}\mathbf{U})\mathcal{D}\bar{\Psi}_{n}\mathcal{D}^{\dagger}$$

$$S_{21} = \mathcal{D}\Psi_{n}\mathcal{D}^{\dagger}\operatorname{diag}(\mathcal{D}\bar{\Psi}_{n}\mathcal{D}^{\dagger}\bar{\mathbf{U}})\operatorname{diag}(\mathcal{D}\bar{\Psi}_{n}\mathcal{D}^{\dagger}\bar{\mathbf{U}})\mathcal{D}^{\dagger}\Psi_{n}\mathcal{D}$$

$$S_{22} = \mathcal{D}\Psi_{n}\mathcal{D}^{\dagger}\operatorname{diag}(\mathcal{D}\bar{\Psi}_{n}\mathcal{D}^{\dagger}\bar{\mathbf{U}})\operatorname{diag}(\mathcal{D}^{\dagger}\Psi_{n}\mathcal{D}\mathbf{U})\mathcal{D}\bar{\Psi}_{n}\mathcal{D}^{\dagger}$$

## Reference

1. K. Kreutz-Delgado, "The Complex Gradient Operator and the CR-Calculus," (2009).

2. J. Nocedal and S. Wright J, *Numerical Optimization*, Springer Series in Operations Research and Financial Engineering (Springer New York, 2006).