Ghost spintronic THz-emitter-array microscope

Terahertz (THz) waves show great potential in nondestructive testing, biodetection and cancer imaging. Despite recent progress in THz wave near-field probes/apertures enabling raster scanning of an object’s surface, an efficient, nonscanning, noninvasive, deep subdiffraction imaging technique remains challenging. Here, we demonstrate THz near-field microscopy using a reconfigurable spintronic THz emitter array (STEA) based on the computational ghost imaging principle. By illuminating an object with the reconfigurable STEA followed by computing the correlation, we can reconstruct an image of the object with deep subdiffraction resolution. By applying an external magnetic field, in-line polarization rotation of the THz wave is realized, making the fused image contrast polarization-free. Time-of-flight (TOF) measurements of coherent THz pulses further enable objects at different distances or depths to be resolved. The demonstrated ghost spintronic THz-emitter-array microscope (GHOSTEAM) is a radically novel imaging tool for THz near-field imaging, opening paradigm-shifting opportunities for nonintrusive label-free bioimaging in a broadband frequency range from 0.1 to 30 THz (namely, 3.3–1000 cm−1).


Section 1. Relevant factors of near-field imaging
Three factors impact the results of near-field ghost imaging, namely sensing distance which determines the spatial resolution, polarization which is relevant to the spatial information and the geometry of DMD which is relevant to the accuracy of coding profile on STEA.

Sensing distance
As sensing distance increasing from near field (≪ λ) to far field, the evanescent waves decay exponentially which results in the larger spatial resolution. In our case, i.e., wavelength λ0 = 600 μm, refractive index n = 1.97 (SiO2), and propagation distance z = 150 nm (thickness of the SiO2 protective layer of our spintronic THz emitter), the critical spatial resolution is estimated as ~188 nm, as described in details below.
We assume that electromagnetic wave with wavelength of λ0 transmits through two slits whose widths are both w and central positions are w apart from each other. The incident electric field is homogeneous with amplitude of 1. Then the field distribution is (1.1) where rect[(x − a)/b] is the rectangle function with central position at a and width of b. The angular spectrum of E 0 (x) can be acquired by Fourier transformation (1.2) where f x represents the spatial frequency. According to the Helmholtz equation, the angular spectrum at z is given by Where k = 2πn/λ 0 = 2π/λ and u x ≡ f x /λ. The electric field at z can be acquired by inverse Fourier transformation (1.4) If the two slits are well resolved, one can find a peak (with peak value of Epeak) in either slit and a valley (with valley value of Evalley) between the two slits. The contrast of |E z (x)| is defined as contrast ≡ (Epeak − Evalley)/Epeak.
(1.5) To quantify the spatial resolution at given propagation distance z (150 nm), the contrast values at slit width varied from 2 nm to 1 μm were calculated using equations (1.1)-(1.5), as shown in Fig. S1a. According to the Rayleigh (Sparrow) criterion, two slits are barely resolved with a contrast equals to 20% (0%), and the corresponding critical slit width was calculated as wR = 188 nm (wS = 100 nm). The critical electric field amplitude |E z (x)| under Rayleigh (Sparrow) criterion is shown in Fig. S1b. Corresponding results with z = 10 nm (blue dashed curves) and z = 1 μm (orange dashed curves) are also given in Figs. S1a and S1b for comparison.
But note that the spatial resolution in our system is ultimately limited by the imaging accuracy of fs-laser (800 nm) on STEA, which is estimated at the scale of micrometers due to the diffraction limit of the 800-nm laser.

Polarization
Polarization has a significant impact on subwavelength imaging. As the feature (slit) size is longer than the wavelength, the transmittance (normalized by slit area A slit ) T(A in /A slit ) of both p-polarized and s-polarized waves approach 1 with decreasing oscillation amplitude, as plotted in Fig. S2a. In Fig. S2b. We calculated the T(A in /A slit ) for polychromatic ppolarized and s-polarized waves (0.1-1.5 THz) within slit width smaller than 600 μm. In subwavelength region (slit width < 60 μm), T(A in /A slit ) of p-polarized wave is increasing sharply while T(A in /A slit ) of s-polarized wave is decaying dramatically, which are consistent with the experimental results (e.g., the region between letters "A" and "E" in Figs. 2b and 2c). Note that the averaged T(A in /A slit ) of two polarizations approaches 1 at most feature size, that is the reason that fused image (Fig. 3d) is more homogeneous.

Impact of diamond micromirror array on coding
Due to the diamond-arrayed micromirrors of the DMD 1 , the accuracy of coding profile on STEA would decrease when a mask pixel is formed with less micromirrors, as illustrated in Fig. S3. It is a fatal factor to our ghost imaging, notably in the case of projected pixel size smaller than 6.5 μm (pixel number corresponds to 128 × 128 in FOV1), as shown in  Wave front is tilted with angle of ~12° in every reflection by on-state DMD and two lenses (focal length of 100 mm) were used to correct the wave front, as described in the main text.
The wave fronts are colored in Rainbow to distinguish their directions. c, One-to-one Image of DMD2 (with some certain mask displayed) at the position of DMD1 acquired by CCD. d,

Section 3. Wave-front measuring and impact on fixed-time ghost imaging
It is critical to make the wave front of pump pulse as flat as possible to cancel out the temporal smearing added by DMDs, since the ghost imaging was operated at some certain fixed time delay.

Wave-front measuring
One can see Supplementary GIFs. 1 and 2 for full raw data and the spatiotemporal THz waveform maps reconstructed with sub-sampled data.

Impact of wave-front tilt on THz ghost imaging
In the premise of guaranteeing the conjugate positions of two DMDs, the optimal correction of wave front is limited by the micromirror tilt angle tolerance of the pair of chosen DMDs, which plays a fatal role in our NGI system. DMD consists of millions micromirrors and each micromirror is switchable between two discrete positions: -α + β and α + β, where α (12°) represents micromirror tilt angle and β (|β| ≤ ~1°) represents the tilt angle tolerance 1 . Both α and β are measured relative to the plane formed by the overall micromirror array (0° flat reference when the micromirrors are parked in their inactive state), as illustrated in Fig. S6a. By comparing the zero-order diffraction angle in the case of normal incidence, the tilt angle differences between every two DMDs among all three DMDs we have were calculated and are shown in the histogram Mask value where Δt represents the time shift acquired by linear fit to peak time delays, c represents the light speed in air and l DMD represents the DMD coding size in the measuring direction (5535 μm in our cases). For a visual comparison, the experimental results whose wave front was corrected by DMD3 and DMD2 are shown in Figs. S6c-S6e. In this case, the time shift in horizontal was calculated as 385 fs × 2 (corresponding deduced tilt angle difference was 2.39° , which is in agreement with the measurement of 2.45° by trigonometry, see Fig. S6b). In horizontal direction (x), the THz transmission decays too fast from the center to the ends to resolve the spatial information of the object, compared to THz transmission in vertical direction (y), as shown in Fig. S6e. It is worth noting that although the time shift can be completely offset by adjusting the position of DMD slightly away from the conjugate position, the mask patterns cannot be projected clearly simultaneously.

Estimation on SNR of ghost image
In practical ghost measurements, the signal vector Y consisting of N elements can be described as = + , (4.1) where O represents the N-length object vector, Φ represents an N × N measurement matrix and e represents an N-length error vector whose each element is associated with each measurement.
The ghost image can be reconstructed by = −1 = − ( + ) = + −1 . It is worth noting that σ e is dependent on the multiplexing method Φ. In our pulsed imaging system, the noises consist of dark noises σ d from detector and source fluctuation σ s among pulses. Here we define two dimensionless parameters γ d ≡ d / 0 and γ s ≡ s / 0 to describe the dark-noise-to-peak ratio and source-fluctuation-to-peak ratio, respectively. In the case of Hadamard multiplexing, half of the incident light is used for each measurement (apart from the pair of first positive and negative measurement), and then σ e,H can be calculated by e, = √ 2 √ 2 + ( 2 )

)
where k represents the number of measurements for each mask and the product term √ indicates that each measurement value in vector Y is a subtraction of one negative measurement from its corresponding positive measurement. While in the case of raster scanning, only E /N of the incident light is used for each measurement (N is usually greater than 1000), so the source fluctuation is negligible compared to the dark noise, and then σ e,R can be calculated by where the number of measurements for each scanning mask is k, for fair comparison in the term of acquisition time. With Equations (4.6) -(4.9), the SNR H and the SNR R can be calculated by , (4.1 ) Here, we simulated the noised ghost images multiplexed using Hadamard matrix with experimental parameters (γ d = 1 −3 , γ s = 7 × 1 −3 , and k = 1 ), which are shown in Fig.  S7. The raster scanning results are also given for comparison. The SNRs calculated from the reconstructed images meet agreement with the corresponding SNRs estimated via Equations (4.10) and (4.11), which are shown at the bottom of Fig. S7b-S7g.

Estimation on the potential frame rate of ghost image
The imaging speed of a GHOSTEAM driven by a MHz oscillator can be enhanced significantly, due to three facts: 1) MHz oscillator can provide minimal acquisition time for single mask (tmask). The shorter tmask is, the faster frame rate (FPS) can be achieved. The minimal tmask is limited by the pulse period t0. The pulse period of oscillator is ~10 ns (e.g., 12.5 ns for an 80-MHz oscillator, see Fig. S8), which is five orders of magnitude smaller than the kHz-amplifier pulse period. But it is worth noting that besides pulse period t0, minimal tmask is limited by the DMD switching time as well, which has typical values of 44-90 μs (depends on the specific DMD type). So, the MHz oscillator can provide minimal tmask at the scale of several hundreds of microseconds. 2) MHz oscillator can provide much more stable illumination for ghost imaging within the identical acquisition time, compared to kHz amplifier. A stable illuminating source can provide ghost images with the same quality (SNR) in much shorter acquisition time. The fluctuation ratio among pulses ("integration time" of 12.5 ns for each data) from oscillator was measured as γs0 = 3.8 × 10 -3 , as shown in Fig. S8. 3) MHz oscillator can efficiently excite the STE to radiate THz pulses with a demonstrated dynamic range of more than 60 dB (1000:1) 3 . The intense illuminating source can also provide ghost images with the same quality (SNR) in much shorter acquisition time.

Figure S8|
Laser pulses from 80 MHz oscillator. The fluctuation ratio among pulses γs0 was measured as 3.8 × 10 -3 from the recorded 400 pulses.
The acquisition time for one ghost image via N-order Hadamard multiplexing equals 2Ntmaskrc, where rc represents the compressive ratio. The frame rate then can be calculated by: Within the "integration time" of tmask, there are γ d ( mask ) = d0 √ 0 / mask , (4.13) γ s ( mask ) = s0 √ 0 / mask , (4.14) Where γ d0 and γ s0 denote the ratio of dark noise to single THz peak and pulse fluctuation ratio within the "integration time" of pulse period t0, respectively. In combination with Equation (4.10), the signal-to-noise ratio of ghost image via Hadamard multiplexing is: The imaging time can be calculated as: GI = mask .(4.17)In the GHOSTEAM driven by a typical 80-MHz oscillator, the pulse period t0 equals 12.5 ns, and the pulse fluctuation ratio γ s0 was measured as 3.8 × 10 -3 , as shown in Fig. S8.
In Ref. [3], the researchers use a 75-MHz oscillator to drive a spintronic THz emitter (pump energy of 4.7 nJ, pump spot diameter of ~10 μm) and obtained THz signals whose spectral dynamic range is well above 60 dB (1000:1) with a maximum around 1.5 THz. The THz signals were detected by a dipole GaAs photoconductive antenna with integration time of 100 ms and the corresponding spectra were not corrected for the detector response. Thus, the dark-noise-to-peak ratio of (1000) -1 can be reached with typical integration time of 100 ms, namely d (1 ) = 1 × 1 −3 .

Boltzmann sigmoidal function
The Boltzmann sigmoidal function subjected to is used to fit the experimental data in Fig. 2d. Three sample curves with different dx (E 1 = , E = 1 and x = ) are given in Fig. S10a. The edge sharpens with reducing dx. According to a 10%-90% criterion of maximum electrical field amplitude 4 , the spatial resolution can be quantified.

Fit to experimental data
In Fig. 2e of the main text, the y-dependent field amplitudes E(y) were individually averaged from the pixels in each row [E(x a ,y), E(x a +Δx,y), …, E(x b ,y)], where x a and x b represent the starting and end pixel position relevant to the slit region, respectively. And the solid curves are reconstructed using piecewise Boltzmann sigmoidal functions with parameters E 1 = 1.9 ( .1 ), E = .1 (1.9 ), dx = . 7 μm and x being a vector consisted of the positions decided by the metal edges. These parameters were derived from the experimental fits to the first four edges, as described below.
According to a 1 %-9 % criterion of maximum electrical field amplitude, the edge length is calculated as ~1 μm when dx equals 0.2276 μm.

Section 6. Data acquisition and processing for THz near-field ghost imaging
Walsh-Hadamard matrix Φ consisting of -1 and 1 was realized by equation To complete a 64 × 64 ghost image, 4096-order Walsh-Hadamard matrix was used to code the masks that displayed on DMD1. The signal vector Y with length of 4096 was measured as where Y p and Y n represent a 4096-lengthed vector consisting of positive mask values and negative mask values, respectively, as illustrated in Fig. S12a. The ghost image can be eventually reconstructed by To minimize the low-frequency amplitude oscillations arise from slow shift of laser power, the negative masks were performed immediately after corresponding positive masks (Fig.  S12b). To suppress the high-frequency amplitude oscillations arise from detector noises and fluctuation among laser pulses, each mask value was averaged by 15 valid data (grey symbols but the ones in rising edge and falling edge, in Fig. S12b). These data were recorded by a lock-in amplifier with integration time of 100 ms. The energy of Y is mainly centralized to mask patterns with low spatial frequencies (e.g., the neighboring masks of mask #1, #128, #256, …, as shown in Fig. S12c).
One can see Supplementary GIFs. 3-5 for the ghost images reconstructed with subsample data. Pulse fluence (repetition rate of 1 kHz) on spintronic THz emitter was 2.88 mJ cm -2 . Each mask was displayed on DMD1 for τDMD = 2 s and each datum (transmitted THz field) was recorded with LIA integration time of τLIA = 100 ms. c, Reference masks of #64, #128, #192, #256 and #320 coded by 4096-order Walsh-Hadamard matrix.

Section 7. Details of image fusion
The fused images are assessed in terms of clearness and directional spatial information amount, respectively.

Clearness of image
Signal-to-noise ratio (SNR) is a widely used assessment parameter for image's quality. The object is known to be binary, and so we choose two regions, one representing a metal area and the other representing a transparent area to calculate the SNR SNR = / , (7.1) where μ is the mean pixel intensity in the transparent region and σ is the standard deviation of the pixel intensity in the metal region. In our experiment, the regions were defined as indicated in Fig. S13a, and corresponding SNR of x-polarized (Fig. 2b), y-polarized ( where X represents the image with L 1 × L pixels, ∇ H and ∇ V are the discretized gradient operators along the horizontal and vertical directions, respectively.

Amount of image's spatial information
The pixels' intensity variation in spatial domain represents the amount of image's spatial information. Due to fact that the THz images were acquired with either T z = | 0 |̂ or T z = | 0 |̂ , there should be a significant difference between horizontal and vertical spatial information amount. Two spatial gradients with mutually orthogonal directions 5 are used to quantify the image's spatial information