Time-Resolved Nonlinear Ghost Imaging

Terahertz (THz) spectroscopy systems are widely employed to retrieve the chemical and material composition of a sample. This is single-handedly the most important driving motivation in the field and has largely contributed to shaping THz science as an independent subject. The limited availability of high-resolution imaging devices, however, still represents a major technological challenge in this promising field of research. In this theoretical work, we tackle this challenge by developing a novel nonlinear ghost imaging approach that conceptually outperforms established single-pixel imaging protocols at inaccessible wavelengths. Our methodology combines nonlinear THz generation with time-resolved field measurements, as enabled by state-of-theart time domain spectroscopy techniques. As an ideal application target, we consider hyperspectral THz imaging of semitransparent samples with nonnegligible delay contribution, and we demonstrate how time-resolved, full-wave acquisition enables accurate spatiotemporal reconstruction of complex inhomogeneous samples.

O ne of the major challenges in photonics is the full-wave characterization of complex field distributions. 1 The possibility of determining the amplitude and phase of the electromagnetic field, in fact, is highly desirable in applications, spanning from biological imaging to material characterization, telecommunications, and optical time reversal. 2−4 As such, it comes as no surprise that imaging devices operating at terahertz (THz) frequencies, where amplitude and phase are experimentally accessible, have been recently subject to intense research. 5 By implementing THz time domain spectroscopy (TDS), in fact, it is possible to perform a direct measurement of the electric field and to quantitatively estimate the refractive index of a sample. 6 Several materials, moreover, possess welldefined spectral signatures in the THz band, which opens the possibility of precisely discriminating their chemical composition. 7 These unique capabilities have potentially disruptive implications in several fields, including deep-tissue biological imaging, security, and single-molecule spectroscopy. 8−10 Despite the large body of potential applications, however, the limited availability of field-sensitive imaging devices remains a major challenge. To date, thermal cameras are routinely employed to characterize THz field distributions. These devices, essentially the equivalent of optical cameras, employ arrays of pyroelectric or bolometric sensors to detect the average intensity distribution of the signal. As such, they are incapable of collecting any information carried by the optical phase, which constitutes the core strength of timeresolved THz detection. To address this limitation, researchers have recently proposed two-dimensional full-wave detectors, composed of arrays of photoconductive antennas or Shack− Hartmann sensors. 11,12 These implementations, however, require complex experimental configurations, and their practicality is still a matter of debate. A valuable alternative to conventional devices can be found in single-pixel imaging, where the sample is illuminated by a well-defined sequence of optical patterns and the resulting scattered signal is collected by a bucket detector. Among different alternatives, ghost imaging (GI) is a very popular and simple single-pixel implementation. 13,14 In its original definition, the incident patterns were randomly generated by a physical process, e.g., through random speckle distributions. In later implementations, however, the computational ghost imaging technique generically refers to a specific imaging process in which the correlated pattern of each realization (i.e., a deterministic pattern) is determined a priori computationally and generated by a proper experimental method. 15−18 Following this approach, it is possible to unambiguously reconstruct the sample by correlating the scattered signal with the spatial distribution of the incident patterns. 19 Remarkably GI has some inherent advantages such as enabling, to some extent, access to resolutions below the diffraction limit or a relative immunity to noisy imaging channels. 20−22 At optical frequencies, single-pixel imaging is usually implemented by measuring the scattered power from the sample. By moving to the THz band, conversely, it has recently become conceivable to combine GI and TDS acquisition, leading to the implementation of a time-resolved, field-based GI protocol where temporal field information can directly contribute to the image reconstruction. 23 Indeed, this combination would lead to a quite different system, in most scenarios only formally related to the original GI. As an immediate example, correlating fields in space−time can be a quite different process when compared to correlating intensity patterns.
In this context, the possibility of performing time-resolved, full-wave analysis of the sample transmission properties represents the desired requirement to characterize a wide range of complex samples, including biological matter, complex three-dimensional scatterers, and semitransparent materials. 24 These systems, in fact, exhibit complex transmission properties that cannot be easily interpreted using a time-independent, intensity-based analysis, which constitutes the standard methodology available today. 25 Despite the large body of research in this field, however, one of the practical challenges in implementing field-sensitive THz GI lies in the realization of controllable structured light. 26,27 Spatial-light-modulator (SLM) devices, which are readily accessible at optical frequencies, cannot be easily extended to the THz domain. Besides, far-field manipulation appears in all cases unsuitable for high-resolution THz imaging, in light of the very coarse diffraction limit stemming from the relatively long THz wavelength. To address these limitations, researchers have implemented indirect patterning schemes (essentially masking devices), such as metallic masks or photoexcited silicon substrates. 25,28,29 While these approaches have been successfully demonstrated in recent THz single-pixel experiments, they are still affected by several conceptual and practical drawbacks, in particular when considering the amplitude and phase properties of the patterned THz beam.
In this work, we devise an entirely new approach to timeresolved, coherent single-pixel imaging. Our methodology, which outperforms standard GI implementations in relevant THz-imaging scenarios, is based on the direct generation of THz structured light through the nonlinear conversion of an optical pattern. Our modeling of the imaging process, which involves an inverse propagation of the temporal traces obtained through TDS detection, demonstrates that nonlinear GI can reconstruct without ambiguities the full response of the sample, including phase and temporal delay contributions. This allows us to obtain a large amount of information on the light−matter interaction between THz illumination and sample. As a result, nonlinear GI represents a complete hyperspectral imaging approach that allows the characterization not only of the sample morphology (with deep subwavelength resolution) but also of its local spectral properties. Our results widen the potential scope of application of time-resolved GI in any part of the spectrum accessible only through nonlinear conversion. These include, among others, thermal and ultrasound imaging, matter waves, and quantum information. 30−32 In the specific framework of THz imaging, our methodology enables the investigation of complex inhomogeneous targets through an innovative single-pixel hyperspectral imaging technique, with significant implications in several key fields and applications, such as deep-tissue biological imaging and time reversal of optical waves. 33,34 ■ RESULTS AND DISCUSSION Time-Resolved Nonlinear GI: Problem Definition. Our main objective is to devise a field-sensitive, single-pixel imaging protocol based on time-resolved measurements of the transmitted electromagnetic field (amplitude and phase).
While inspired by standard GI (which relies on timeindependent intensity analysis), our methodology is enabled by the assumption that the field distribution is determined by the optical intensity distribution impinging onto a THznonlinear converting element. In this sense, this scenario fits the standard requirements of high-energy THz TDS systems exploiting thin nonlinear crystals.
We model the transmission properties of the sample by defining a spatiotemporal transfer function, T(x, y, t). The explicit dependence on time, which is pivotal in the accurate reconstruction of a sample with complex transmission properties, allows us to retain the full amount of temporal effects impressed by the object on the incident pattern. Without loss of generality, we express the electromagnetic field impinging on the object as E n − (x, y, t) = P n (x, y) f(t) where P(x, y) and f(t) represent the spatial pattern and the temporal profile of the THz pulse and where we assumed that the incident pattern is not affected by spatiotemporal coupling. Space-time separability is a standard requirement in computational imaging, which assumes a precise knowledge of the spatial information encoded in the incident pattern. We stress that if the THz pattern originates from the conversion of an optical pattern, this condition is realized in realistic experimental scenarios. In terms of the electromagnetic field, the transfer function of the object T(x, y, t) is reconstructed as follows: n n n n n n n (1) where ⟨···⟩ n represents an expected value over the pattern ensemble and C n (t) are expansion coefficients defined as (2) * being the convolution operator. In analogy with standard GI implementations, the coefficients C n (t) coincide with the integrated scattered field measured by the bucket detector. 35 It is worth stressing, however, that eq 2 involves a direct measurement of the electromagnetic field E n + (x, y, t) = T(x, y, t)*E n − (x, y, t), in stark contrast with the intensity-based measurements of standard GI. A direct measurement of the scattered field, which is rather inaccessible in the optical domain, can be performed at THz frequencies by employing a TDS detection scheme. Moreover, eq 2 provides significant advantages when designing actual experiments. The C n defined in eq 2, in fact, coincide with the k x = k y = 0 spectral components of the THz transmitted beam (Fraunhofer approximation).
One of the main challenges when translating eqs 1 and 2 in real experiments, however, lies in the generation and control of the incident field distribution E n − (x, y, t). This is particularly true at frequencies where standard beam-shaping equipment is not easily accessible. We tackled this problem by considering a structured wave generated through a nonlinear optical process, as established, for example, in the state-of-the-art generation of THz radiation from ultrafast optical illumination. 36,37 As the THz beam is generated through a quadratic transformation (difference-frequency generation), there is a local linear dependence between the THz field and the incident optical intensity, namely, E THz (x, y, t) ∝ I opt (x, y, t). Such a relation allows precise control over the THz field profile both in time and space by shaping the incident optical pulse. With this approach, any optical pattern generated through a standard SLM device can act as a direct source of THz structured

ACS Photonics
Article radiation, with the advantage that the spatial resolution of the THz patterns can be pushed way below the standard THz length scale, as it is only bound by the diffraction limit of the optical beam. Even in these conditions, the individual pixel intensity and the phase profile of the patterned THz beam remain completely independent from its spatial distribution, and they depend only on the amplitude and temporal delay of the incident optical beam. Indeed, this is a requirement for the implementation of time-resolved, field-based GI methodology, and it is very difficult to achieve through metallic masks or optically induced photocarrier masks. 25,38 Masks with highresolution features, in fact, are bound to the stringent physics of subwavelength apertures, where the transmitted amplitude and phase depend directly on the aperture size. 39 In these implementations, the relationship between the incident optical pattern and the structured THz beam involves the Green's function of the subwavelength apertures composing the mask. In our nonlinear approach, conversely, the relationship between incident optical intensity and THz beam is of direct proportionality, with the additional advantage that the optical beam does not experience any spatiotemporal coupling on the surface of the generating crystal. Therefore, nonlinear frequency conversion represents a suitable method to generate subwavelength THz patterns characterized by space−time separability, a critical requirement of computational imaging. Moreover, this approach allows us to generate THz patterns whose intensity scales as 1/r 2 (as opposed to the well-known 1/r 6 scaling of subwavelength apertures). 40−43 As a result, nonlinear frequency conversion of deeply subwavelength optical patterns allows generating THz structured light without compromising its spatial resolution, which is limited only by the "focusing power" of the fundamental beam (up to ∼ λ λ 2 600 optical THz ). 44,45 Besides these conceptual challenges, opti-cally induced transient photocarrier masks, which represent the most developed benchmark in the field of THz GI, require careful consideration of additional elements. For example, the masking effect induced by transient photocarriers is not restricted to the semiconductor surface in light of the relatively large skin depth at THz frequencies (e.g., several μm in photoexcited silicon 38 ). The excited pixel in the mask always represents a dark THz pixel, with any residual THz transmitted reducing the pattern contrast. This makes this implementation agile mostly with static binary patterns. By comparison, a nonlinear generation scheme potentially allows for analogic control of the generated THz profile (i.e., "gray" scales), in terms of both THz amplitude and temporal delay. Reference Nonlinear GI Experimental Setup. Our reference nonlinear GI embodiment is illustrated in Figure 1. A structured optical beam E opt (x, y, t) (red) is generated through a standard SLM device operating at optical frequencies ( Figure 1a). Upon impinging on the surface of a generating crystal slab of thickness z 0 (e.g., ZnTe), the optical pattern is nonlinearly converted to THz frequencies, generating a structured THz field, which we denote as E THz (x, y, t) (purple). Upon generation, the THz pulse propagates across the crystal and interacts with the sample (Figure 1b), which we assume for simplicity to be located at z = z 0 . After the interaction, the spatial spectral component k x = k y = 0 (the spatial mode propagating collinearly with z) of the transmitted THz field is collected through a single-pixel detector placed on the Fourier plane and analyzed by a standard TDS setup. In practical terms, these components can be measured by placing the object in the focal point of a lens and by acquiring the signal in the central point of the opposite focal plane. The possibility of implementing a point-like detection in the Fourier plane has critical consequences on the imaging capabilities of our setup. Standard GI implementa-

ACS Photonics
Article tions, in fact, require that the total scattered is collected and integrated by the bucket detector, and, as a result, the numerical aperture of the latter has important consequences on the maximum achievable resolution. Roughly speaking, the bucket detector must be able to acquire as much scattered field as possible to be accurate. By performing a point-like detection, this specific aspect is not directly relevant in the reconstruction approach, although in this case the final signal-to-noise-ratio (SNR) of the measurements is determined by the spatial spectrum of the patterns. A crucial consideration when implementing single-pixel imaging schemes lies in the choice of the incident patterns. In our analysis, we implemented an imaging protocol based on the Walsh−Hadamard encoding, which is known to maximize the SNR. 46 The Walsh patterns form an orthogonal set of binary matrices, which can be employed to expand any two-dimensional images in a deterministic fashion, opening to the implementation of advanced compressed sensing algorithms based on the sample sparsity. 47 Spatiotemporal Coupling in Deeply Subwavelength THz Structured Light. When dealing with subwavelength THz patterns, one of the major practical challenges is posed by the distance between the source plane and the object, as thoroughly discussed in ref 38 and references therein. Such a constraint usually limits the effective resolution achievable by the system, as subwavelength THz patterns experience diffraction while propagating inside the substrate due to spatiotemporal coupling. To quantitatively illustrate this point, we express the structured optical field E ⃗ n opt (x, y, t) as follows: where H n (x, y) is a Walsh−Hadamard pattern, f(t) is the temporal profile of the pulse, and p̂is the beam polarization unit vector, which allows us to retain any polarization effect induced by spatiotemporal coupling at subwavelength propagation distances. On the generation plane (the crystal surface), the optical beam induces a nonlinear polarization field of the form where χ (2) is the nonlinear susceptibility of the crystal and where we assumed that the polarization field is copolarized with the incident beam to retain a lighter notation (the extension to a fully vector χ (2) tensor is straightforward). The nonlinear polarization from eq 4 generates a THz field defined as follows: where N is a normalization constant, which we will assume, without loss of generality, to equal unity. The generated THz beam propagates within the material according to THz THz (6) where G ⃡ (x, y, z, t) is the dyadic Green's tensor (see Methods). The use of the dyadic Green's tensor, as opposed to the standard Green's function of scalar diffraction theory, is a strict requirement to devise a vector model of spatiotemporal coupling upon propagation. We describe the transmitted field due to the interaction between the THz pulse and object as follows:

ACS Photonics
n n 0 0 where E ⃗ n − is the field immediately before the sample, E ⃗ n + is the field immediately after, and ϵ is an infinitesimal spatial displacement (namely, ϵ → 0). With this position, the experimental TDS signal of the transmitted field simply reads: where p TDS is the detection polarization state, which can be controlled, for example, by rotating the ZnTe crystal in a standard electro-optic detection. By substituting eqs 7 and 8 into eq 2, the single-pixel amplitudes C n (t), coinciding with the TDS signal measured on the Fourier plane, are rewritten as follows: n n TDS THz (9) As can be evinced by comparing eqs 2 and 9, the propagation within the crystal leads to a rather different physical scenario. In the presence of subwavelength illumination, in fact, the object does not directly interact with the incident pattern, but with a propagated version of it. Diffraction and spatiotemporal coupling, in fact, can ultimately lead to a significant corruption of the initial pattern, as illustrated in Figure 2, which showcases the propagation of a subwavelength Walsh pattern (generated by nonlinear frequency conversion) as a function of the propagation distance z. While the detrimental effect of spatiotemporal coupling could be mitigated by reducing the distance between the generation plane and the object, this is a significant practical constraint. In addition the consequent reduction of source thickness in general leads to a strong decrease in nonlinear conversion efficiency and represents a technological challenge. 48 In practical subwavelength imaging implementations, there is always a relevant propagation distance between source plane and object, especially if an interaction between the optical field and the object is undesirable, which implies an optical depleting or reflecting region. For convenience we will assume that this situation is realized in the nonlinear generation crystal. We argue that the access to the full temporal evolution of the transmitted field is a requirement in time-domain hyperspectral imaging and enables the reconstruction of the sample even in the presence of thick substrates. Following eqs 7−9, the reconstructed transfer function in the presence of a thick substrate reads as follows: n n n n n n n (10) As such, the reconstructed transfer function T′ involves not only the transmission properties of the object but also the propagation across the substrate. Formally, this corresponds to a reconstructed transfer function of the form T′ = T*G. An interesting question is whether the latter can be numerically inverted to extract the sample response.
Time-Resolved Image Reconstruction in a Noisy Environment. To overcome any limitation posed by spatiotemporal coupling through the generation crystal, we implemented an image reconstruction scheme that takes into account the full amount of temporal information provided by the TDS detection. This leads to quite different outcomes when compared to the state-of-the-art THz GI, where the imaging operates considering fixed TDS delays separately ("fixed-time" GI). In the presence of strong spatiotemporal

ACS Photonics
Article coupling, as in the case of deeply subwavelength structured light, the knowledge of the field at a fixed time is very limiting in extracting the sample morphology. In addition, the temporal response reconstructed simply applying the reconstruction at different delays cannot be used to completely determine the optical properties of the sample, in particular when considering objects with a complex spectral response or composition. To illustrate this point, we integrated eqs 3−9 in the Fourier domain. Our numerical results are shown in Figure 3. As a relevant case study, we characterized the transmission from a two-dimensional gold mask. The system was probed through a full set of Walsh pattern (Hadamard order N = 32), propagating across a crystal of thickness z 0 = 117 μm. To quantitatively assess the reconstruction capabilities of our approach, we measured the correlation between reconstructed and original samples by implementing a structural similarity index method (SSIM) analysis. 49 The SSIM was recently proposed as an extension to the established Pearson correlation coefficient, and it provides a detailed comparison in terms of luminosity (mean value), contrast (variance), and structure (covariance) of two spatial distributions (see Methods section for further details). To assess the robustness of our approach to experimental noise, we also included a white-noise term in the TDS trace with a statistical SNR of 5% and compatible with state-of-the-art experimental conditions. In Figure 3a−c we first considered an imaging target with spatial resolution larger than the incident wavelength (λ = 234 μm, target resolution = 468 μm).
In these conditions, the SSIM reaches a maximum value of 0.84 (Figure 3a, orange line) in correspondence with the peak transmission. Intuitively, all the SSIM peaks are placed in correspondence with the maximum field, i.e., in the conditions of minimum SNR. When considering a deeply subwavelength object, however, fixed-time measurements do not provide a suitable morphological reconstruction at any time. This is illustrated in Figure 3d−f, where we considered a target resolution of 20 μm, corresponding to approximately λ/12. Differently from the previous case, the reconstructed image at the SSIM peaks has a rather poor similarity to the original sample (Figure 3e,f).
Even in these conditions, however, it is possible to reconstruct the sample by adopting an entirely different approach: rather than considering the transmitted signal at specific times, we took into account the full amount of information provided by the time-resolved signal. Such information can be processed through a spatiotemporal inverse-propagation operator inspired by the theory of Wiener filters. In the Fourier domain, we define the inversepropagation operator as follows: where G ⃡ (k x , k y , k z , ω) is the dyadic Green's function in the Fourier domain, NSR = 1/SNR is the average spectral noiseto-signal ratio, and α is a fitting parameter (α = 0.1 in our simulations). The action of the operator W ⃡ is to inversepropagate the experimental TDS trace through the crystal and reduce the effects of spatiotemporal coupling (see Figure 4a). In the absence of noise, the operator W ⃡ coincides straightforwardly with the inverse Green's tensor. In the presence of experimental noise, conversely, the NSR term in the denominator of eq 11 is a strict requirement, as in these conditions the transfer function cannot be univocally reconstructed with a standard inverse filter. 50 We verified the validity of our approach by applying the operator from eq 11 to the cases of Figure 2. The reconstructed images are shown in Figure 4b,c. Remarkably, even in the presence of subwavelength illumination ( Figure 4c) the inverse-propagation filter demonstrates a remarkable capability of reconstructing the sample morphology, in sharp contrast with the fixed-time reconstruction (cf. Figure 3d−f). The demonstrated superiority of time-resolved image reconstruction based on inversepropagation holds also in the absence of experimental noise, as illustrated in Supplementary Figure 1. It is worth mentioning that the feasibility of an inverse-propagation transformation strongly relies on the coherence properties of the THz incident pattern generated by nonlinear frequency conversion. When compared to the case of optically induced transient photomasks with subwavelength features, in fact, the THz pattern is defined as the convolution between the incident THz beam and the Green's function of the subwavelength mask. The latter is intrinsically pattern-dependent and hard to model in a general way due to the complex space−time light−matter interaction occurring between the incident beam and the excited photocarriers.
Hyperspectral Single-Pixel Imaging at THz Frequencies. The main properties of our nonlinear time-resolved approach provide significant advantages when considering complex inhomogeneous samples. It is widely acknowledged that one of the main advances enabled by TDS analysis is the

ACS Photonics
Article capability to accurately reveal the spectral response of a sample. 7,23 An intriguing question is whether time-resolved nonlinear reconstruction could be employed to perform hyperspectral imaging, i.e., to access the local spectral information for each spatial point of the reconstructed image. Such kind of measurement, if feasible, would be clearly pivotal in spectral analysis and material recognition of complex semitransparent samples, which currently pose a great challenge to established single-pixel protocols. We verified this possibility by considering a semitransparent target of thickness δ and composed of three different materials: gold, Kapton, and Teflon ( Figure 5, inset). Kapton and Teflon exhibit a rather different spectral response at THz frequencies (in particular in the 1−5 THz region), where Teflon is a THztransparent polymer and Kapton is strongly absorbing. As such, the choice of this materials allows for a straightforward discrimination of the different materials composing the object and their spatial distribution.
In terms of its material composition, the transfer function for the object can be expressed as follows: where τ FP is a transfer function which depends on the local refractive index distribution n(x, y, t) = n(x, y, t) + ik(x, y, t).
Note that in our analysis we took into account the full material dispersion of the target. 51 Without loss of generality, we restricted our analysis to the case of a thin object with features comparable with the THz wavelength and with a thickness of δ = 10 μm. Under these conditions, the transmission function τ FP can be easily expressed in terms of a Fabry−Peŕot response (see Methods section for further details).
The results of our analysis are shown in Figure 5. As illustrated in Figure 5a−c, the hyperspectral image obtained by inverse-propagating the reconstructed signal allows easy discrimination of the different materials in the object. Since the Teflon portion of the sample is quite transparent for frequencies below 5 THz, in this region the hyperspectral image is composed mainly of Teflon components (Figure 5d). Conversely, for frequencies between 6 and 8 THz, where Teflon is mostly opaque and Kapton is transparent, the hyperspectral image corresponds to the Kapton portions of the object. The high degree of chemical selectivity is confirmed by SSIM analysis, whose results are reported in Figure 5e. It is worth mentioning that, despite the drastically different spectral response of the two materials, the transmission spectra from the very thin sample look rather similar (Figure 5d).

■ CONCLUSIONS
In this work, we demonstrated a novel approach to deepsubwavelength single-pixel imaging based on nonlinear pattern generation and time-resolved field measurements. As a relevant case study, we considered single-pixel imaging at THz frequencies, where time-resolved field measurements are readily available. Interestingly, the possibility of performing time-resolved measurements directly descends from our nonlinear pattern generation scheme, which provides an advantageous mechanism to produce structured light at THz frequencies. Time-resolved measurements, moreover, allowed

ACS Photonics
Article us to perform sample reconstruction even in challenging conditions, including distant objects and semitransparent inhomogeneous targets. In particular, we were able to demonstrate the feasibility of single-pixel hyperspectral imaging, which paves the way to a new methodology in twodimensional material characterization.
By comparing our methodology with the state-of-the-art at this frequency range, we have demonstrated the conceptual advantages of a coherent time-resolved measurement in the reconstruction of complex inhomogeneous samples. Interestingly, we have shown how accurate spatial and spectral reconstruction of the sample can be achieved only by inverse propagating the time-resolved TDS traces, i.e., by eliminating the effects of spatiotemporal coupling through the nonlinear generating crystal. Such a result constitutes, in our opinion, a valuable extension of current methodologies and allows the systematic investigation of hyperspectral single-pixel imaging at THz frequencies. Interestingly, the performance of our nonlinear GI protocol could benefit from combining singlepixel field detection with established super-resolution imaging techniques, as in the case of optical superoscillation imaging. 52 We believe that the demonstration of a time-resolved singlepixel methodology based on full-wave measurements has profound implications that go beyond the field of THz imaging. Single-pixel detection and ghost imaging, for example, play a crucial role in the development of advanced applications in the infrared and microwaves. In these spectral regions, nonlinear frequency conversion and bucket detection are actively employed in the realization of quantum radar and sensing technologies. 53 At the same time, it has been recently shown that time-resolved intensity measurements can be employed to achieve quite remarkable performances in deeptissue biological imaging at optical and infrared frequencies, such as single-shot imaging through thick and opaque samples. 3,54 These applications rely on reconstructing the full-wave transmission properties of complex samples, such as turbid media or biological tissue. An intriguing question is whether such methodologies could be implemented at THz frequencies, where amplitude and phase are directly measured.

■ METHODS
Nonlinear Generation and Propagation of THz Subwavelength Patterns. We generated structured THz pulses by nonlinear optical rectification of a structured optical pulse. In our analysis, we considered an optical pulse with a Gaussian profile, centered at λ = 800 nm and with a duration of 300 or 80 fs used in the time-resolved GI and hyperspectral imaging reconstruction, respectively. Upon generation on the crystal surface, the THz structured pulse propagates according to eq 6, where the dyadic Green's function G ⃡ (x, y, z, t) is defined as 55 where ε THz = 13.7 is the dielectric constant of the crystal at THz frequencies (we assumed a value compatible with a THz quadratic nonlinear crystal) and Θ(t) is the Heaviside step function. Structural Similarity Index Method. We quantitatively assessed the correlation between the reconstructed images and the original sample by employing the SSIM, initially proposed in ref 49. Standardly, the spatial correlation between images is measured by the Pearson correlation coefficient ρ P , which is defined as ρ σ σ σ = x y ( , ) xy x y P (15) In the SSIM model, conversely, the similarity between two images is defined in terms of the luminance (l), contrast (c), and structure (s) coefficients, which are defined, respectively, as follows: x y x y x y x y 2 2 2 2 P (16) In terms of these coefficients, the SSIM similarity index is defined as = α β γ l c s SSIM (17) where α, β, and γ are fitting weights. In this analysis, we considered the following set of exponents: α = 0.01, β = 0.01, and γ = 2. Such a choice allows us to underline the structural information provided by the Pearson coefficient. From our analysis, we considered a minimum similarity index SSIM of 14%, corresponding to the SSIM between the original image and a nonstructured Gaussian field distribution.
Transfer Function for Inhomogeneous Samples. We modeled the transmission properties of our semitransparent inhomogeneous sample in terms of a Fabry−Peŕot response. 57 In the frequency domain, this corresponds to an analytic transfer function of the form  (18) where φ(x, y, ω) = n(x, y, ω)ωδ/c, n(x, y, ω) is the local refractive index distribution of the sample, and we defined the transmission τ and reflection ρ coefficients as

ACS Photonics
Article