**4. Performance of Field-Based Ghost-Imaging Detection in the Fourier Plane**

The possibility of performing field-sensitive detection provides TIMING with a significant advantage when compared with traditional GI. However, the typical GI correlation between detection parameters and image fidelity is broken by the nonlinear ghost imaging transformation, i.e., the need for establishing a correlation between coherent-field detection and the optical intensity patterns. More precisely, the implementation of a field average in the image extraction radically changes the way the image quality depends on the experimental parameters. Standard GI reconstruction relies on detecting the integrated scattered field to estimate the spatial correlation between the incident patterns and the sample, where:

$$\mathcal{L}\_n = \int \mathrm{d}\mathbf{x} \mathrm{d}y \, I\_n^+(\mathbf{x}, y) = \int \mathrm{d}\mathbf{x} \mathrm{d}y \mathrm{d}t' \Big| T(\mathbf{x}, y, t - t') E\_n^-(\mathbf{x}, y, t) \Big|^2. \tag{7}$$

This corresponds to the direct acquisition of the total scattered field with a standard bucket detector, which integrates the transmitted intensity distribution. Fundamentally, it is an estimator of the total scattered power, and as such, it is directly affected by the numerical aperture of the detector and by the distance between the detector and the sample. As discussed in the literature on optical GI, both these factors directly fix the amount of information that is available when reconstructing the image and directly affect its fidelity [15].

TIMING inherits the direct detection of the scattered THz field distribution from time-domain spectroscopy systems. By operating directly on the electric field, it allows for measuring the average THz scattered field (in a fully coherent sense) by performing a point-like detection in the Fourier plane. As defined by Equation (3), the coefficients *Cn* can be obtained by measuring the *kx*, *ky* = 0 spectral components of the THz transmitted field:

$$\mathcal{C}\_{n}(t) \;= \int \mathbf{dx} \mathbf{dy} \, E\_{n}^{+}(\mathbf{x}, \mathbf{y}, t) \;= \mathcal{F} \left[ E\_{n}^{+}(\mathbf{x}, \mathbf{y}, t) \right] \; \vert\_{\mathbf{k}\_{x} = 0, k\_{y} = 0} \tag{8}$$

This implementation implies that the experimental measurement of the correlations *Cn* is not limited at all by the numerical aperture of the bucket detector. This type of measurement can be obtained by placing the object in the focal point of an arbitrary lens and by acquiring the signal in the central point of the opposite focal plane (Figure 1a). The electric field in the focal plane reads as follows:

$$E\_{\text{focal}}(\mathbf{x}', \mathbf{y}') \propto \mathcal{F} \left[ E\_n^+(\mathbf{x}, \mathbf{y}, t) \right] \left( k\_{\mathbf{x}} \;= \; \frac{\mathbf{x}'}{\lambda f}, k\_{\mathbf{y}} \;= \; \frac{\mathbf{y}'}{\lambda f} \right) \tag{9}$$

where *x* and *y* are the physical coordinates in the Fourier plane [41].

However, in terms of implementation, the detector samples a finite small area of the Fourier plane with an area-sampling function *PH kx*, *ky* , obtaining the estimation *Cn* (*t*):

$$\mathcal{L}\_{n}^{'(t)} = \int \boldsymbol{P} \boldsymbol{H} \Big(\boldsymbol{k}\_{\mathbf{x}}, \boldsymbol{k}\_{\mathbf{y}}\big) \ast \, \mathcal{F} \Big[\boldsymbol{E}\_{n}^{+}(\mathbf{x}, \mathbf{y}, \mathbf{t})\Big] \, d\mathbf{k}\_{\mathbf{x}} d\mathbf{k}\_{\mathbf{y}}.\tag{10}$$

where *PH kx*, *ky* is physically represented by the profile of the probe beam in the electro-optical sampling (e.g., a Gaussian function), or by the shape of any aperture implemented in front of the nonlinear detection to fix its interaction area with the THz field.

The accuracy of the measurements is then directly related to how "point-like" our detection can be made. Although one could be tempted to foresee a general benefit of the high signal-to-noise ratio (SNR) resulting from large detection apertures as in the standard GI, this is also a source of artefacts, fundamentally establishing a trade-off between SNR and fidelity.

Figure 3 illustrates the effects of the size *d* of the sampling function *PH x* = *kx*λ*f*, *y* = *ky*λ*f* on the image reconstruction fidelity (Figure 3e). Interestingly, the reduction of fidelity observed for increasing the sampling diameter is different from the typical limitations in standard imaging. In our case, a too-large area sampling function in the Fourier plane did not lead to a reduction in the discernible details but rather in the disappearance of entire parts of the image (see Figure 3e, insets).

Similarly, in Figure 4, we illustrate the effect of a misalignment of the sampling function *PH* centre with respect to the centre of the Fourier plane. Trivially, the spatial correlation between the reconstructed and original images peaks at the centre of the Fourier plane and swiftly decayed in the case of off-axis detection (Figure 4a). In these conditions, the reconstructed image showed the appearance of spurious spatial frequencies, corresponding to the *kx*, *ky* sampling position (Figure 4b,d). Interestingly, however, the overall morphology and details of the image were still present in the images, and no noticeable blurring occurred.

**Figure 3.** Influence of the pinhole size on the Fourier detection of TIMING reconstruction coefficients. **(a**–**d**) The spatial average of the transmitted field (**b**) associated with each incident pattern (**a**) could be measured by performing a point-like detection in the centre of the Fourier plane (**c**,**d**). In realistic implementations, the centre of the Fourier plane is sampled using a sampling function *PH* of finite diameter *d*. (**e**) Spatial correlation between the reconstructed and original image as a function of the sampling function diameter. A departure from the point-like approximation led to a significant corruption of the reconstructed image (insets). Interestingly, the typical image degradation did not necessarily involve the total disappearance of highly resolved details.

**Figure 4.** Influence of the pinhole displacement on the Fourier detection of TIMING reconstruction coefficients. (**a**) Spatial correlation between the reconstructed and original image as a function of the sampling function position in the focal plane. The displacement*(*Δ*x,* Δ*y)*was measured with respect to the lens axis and the sampling function diameter was set to *d*=0.36 mm, corresponding to a spatial correlation of 100% at the centre of the Fourier plane (cf. Figure 3e). (**b**–**d**) Examples of image reconstruction with off-axis detection, illustrating the appearance of spurious spatial frequencies. Interestingly, the object morphology was still noticeable, even at a relatively large distance from the optical axis.
