Subwavelength terahertz imaging via virtual superlensing in the radiating near field

Imaging with resolutions much below the wavelength λ – now common in the visible spectrum – remains challenging at lower frequencies, where exponentially decaying evanescent waves are generally measured using a tip or antenna close to an object. Such approaches are often problematic because probes can perturb the near-field itself. Here we show that information encoded in evanescent waves can be probed further than previously thought, by reconstructing truthful images of the near-field through selective amplification of evanescent waves, akin to a virtual superlens that images the near field without perturbing it. We quantify trade-offs between noise and measurement distance, experimentally demonstrating reconstruction of complex images with subwavelength features down to a resolution of λ/7 and amplitude signal-to-noise ratios < 25dB between 0.18–1.5 THz. Our procedure can be implemented with any near-field probe, greatly relaxes experimental requirements for subwavelength imaging at sub-optical frequencies and opens the door to non-invasive near-field scanning.


INTRODUCTION
The diffraction limit is a consequence of the evanescent decay of high spatial frequencies in standard materials [1].Conventional imaging techniques that collect light far from an object are typically bound by this limit, and much effort has been invested in developing ways to overcome it.Many techniques now provide resolutions well below the diffraction limit [2,3], relying either on near-field probing through a scanning tip [4,5], stochastic sets of scatterers or fluorophores in the immediate vicinity of the object to be imaged [6,7], or nonlinear effects [8,9], only really possible in the optical spectrum.Methods to reconstruct sub-diffraction details from linear far fields also exist, but typically require some prior knowledge or assumptions on the nature of the object [10].
Lower frequency sub-wavelength imaging techniques (e.g., GHz, THz) typically rely on scanning antennas in an object's near field [11,12].Imaging at terahertz frequencies (0.1-10 THz) would particularly benefit from any improvement in the ability to image below the diffraction limit [13], due to its many applications in biomedicine [14][15][16][17], which is hindered by established experimental challenges [18].We refer the reader to Ref. [13] for a recent review on recent developments in terahertz imaging techniques.The state-of-the art, in terms of resolution (<100 nm), is given by nanoscale THz scanning probe microscopy [19], a sophisticated technique that can only access small-area planar surfaces.A practical alternative which can be incorporated into any commercial THz time domain spectroscopy setup [20], is provided by near-field photoconductive detector antennas [12,21], which directly probe both amplitude and phase over cm 2 areas with resolutions of order 10 µm at the site of the antenna.However, under standard laboratory conditions, it is common for such antennas to be hundreds of micrometers away from the object to be imaged [22,23].Such large measurement distances are arguably desirable, since local fields can be heavily perturbed by the presence of near-field probes [24].Furthemore, such antennas are expensive and delicate, making near-contact scans a risky proposition.At THz frequencies, such distances are associated with the radiating near field region (i.e., distances between λ/2π and λ from the object [25]), where high spatial frequencies decay significantly, preventing genuine sub-wavelength imaging [1].
Ideally, the exponential decay of the evanescent field would be avoided or reversed.This can be achieved with so-called superlenses (SL) [26,27] and hyperlenses (HL) [28], which respectively amplify or propagate the evanescent fields, and which have been demonstrated over much of the electromagnetic spectrum [11,27,[29][30][31][32][33].Such approaches still carry two challenges: (i) the resolution of the highest spatial frequencies is adversely affected by even modest losses [34]; (ii) most geometries transfer the near field information without spatial magnification which converts evanescent waves into propagating waves [30,35].Fields must still be measured in the near field of the SL or HL [36], shifting the problem rather than solving it.At THz frequencies, anisotropic metamaterials have been used for subwavelength propagation of near-field information across finite slabs [30,37,38] using Fabry-Perot resonance-induced evanescent amplification [39,40], but a THz SL which truly amplifies decayed evanescent fields has so far been beyond reach.
In the absence of noise, reversing evanescent decay does not necessarily require physical devices, but can be done numerically instead.Evanescent waves contribute to the field measured at any distance, superimposing high spatial frequency fluctuations on top of low spatial frequency radiating fields.If the fluctuations can be resolved, they can also be amplified to regain the original field.To our knowledge, such a scheme has not been considered, presumably because high-spatial frequencies commonly decay well below instrument noise.Here we show that practical noise levels allow for a useful extraction of decayed information in the radiating near field region, also providing a way for increasing the resolution in the reactive near field region.In effect, we experimentally demonstrate a virtual superlens through post-processing, reconstructing previously indiscernible sub-wavelength spatial features contained within complex images, with demonstrated resolution down to λ/7.Our approach is general, provided that lownoise, phase-resolved fields can be measured.This opens the possibility of measuring near fields without perturbing them -which would be particularly useful when resolving fields in structures that are sensitive to perturbations, e.g., high-Q/topological resonances [41,42] and photonic crystal defects [43].

Approach and implementation
Figure 1 shows a schematic of our approach, which aims to image a planar source object possessing subwavelength features with a field E obj (x, y, z = 0).The total field at r = (x, y, z) is given by a Fourier expansion [26] where σ sums over polarizations, k = (k x , k y , k z ), Ẽσ can be obtained from the Fourier transform of E obj (x, y, z = 0), and k 0 = 2π/λ is the free space wavenumber.Propagating waves (Fig. 1, dashed) carry information emerging from spatial frequencies satisfying k 2 x + k 2 y < k 2 0 and impose a lower limit on the spatial features d which can be resolved in the far field.Evanescent waves (curves) carry sub-wavelength spatial frequencies satisfying k 2 x + k 2 y > k 2 0 and exponentially decay in free space.As a result, the fine details of an image possessing spatial features d λ, detected by a near field probe at a distance z = L, cannot be resolved.This process can be straightforwardly reversed via the transformation (x, y, z) → (−x, −y, −z) over a subsequent length L, by numerically reversing the phase accumulated by the propagating waves and amplifying the evanescent waves (green curves in Fig. 1).In practice, we measure an x-polarized field E x (x, y, z = L).Starting from Fourier components of the measured field Ẽx (k y , k z ), the electric field after the virtual SL is given by: where k z ∈ C follows Eq. ( 2), with arbitrarily high spatial frequencies.For large spatial frequencies, k z is imaginary and leads to exponential amplification exp |k z |L.This process is equivalent to "superlensing" [26], achieving the same effect (Fig. 1, blue dotted line).Simply reversing the phase without amplifying evanescent waves, that is limiting the integrals to real k z , is akin to what an ideal conventional (far field) lens achieves ("lensing").In Eq. (3) higher spatial frequencies lead to larger amplification terms so that, after amplification, high spatial frequency noise is bound to dominate over any signal at lower spatial frequencies.A simple numerical example showcases the issue: We consider 2D finite element method calculations where the domain is infinite in y, with TM polarized fields (non-zero magnetic field in y). Figure 2(a) shows |E x | emerging from two perfectly conducting apertures with width-and edge-to-edge-separation of d = λ/5.At a distance z = λ/10, the two apertures can be distinguished in the field (blue curve in Fig. 2(a)).At a distance of z = λ/2 (red), however, this is no longer the case.Figure 2(b) shows the associated spatial Fourier transforms.The minimum magnitude of k x required in order to resolve d is given by k x /k 0 = λ/2d = 2.5.While the source's spatial Fourier spectrum extends beyond k x = 10k 0 , at z = λ/2 (i.e., in the radiating near field region) most field components with spatial frequencies inside |k| < k max have exponentially decayed 20 dB below that of propagating waves, so that the apertures cannot be distinguished.Note that while the source spectrum (Fig. 2(b), blue) is smooth, the radiating near field spectrum (Fig. 2(b), red) presents numerical noise starting even at modest k x /k 0 values.Using this noise level compared to the maximum signal, the amplitude signal-to-noise (SNR) is ∼ 30 dB.
We now implement the superlens procedure given by Eq. ( 3) to the complex field associated with the red curves in Fig. 2(a) and 2(b) obtained at z = λ/2.The result is shown as a black dotted line in Fig. 2(c).The image reconstruction of the two apertures has clearly failed: the associated spatial Fourier spectrum (black line in Fig. 2(d)) shows that noise at high spatial frequencies has been amplified to exceed the amplitude of any signal, including that of the propagating waves.However, not all is lost: comparing the black curve amplified spectrum with the blue curve in Fig. 2(b), the amplified spectra match with the original for k x /k 0 2.5, where the signal in the evanescent spectrum was originally above the noise level (Fig. 2(b), purple).We thus apply a spatial low-pass filter in k-space after superlensing, setting all |k x | > 2.5k 0 to zero.The result is shown as blue curves in Fig. 2(c,d), where we find that the virtual SL now resolves the apertures as desired.In contrast, filtering all non-propagating waves as per a conventional lens, i.e., setting regions where k > k 0 to zero (green curves in Fig. 2(c,d)), does not allow us to resolve the apertures, even though the phase is reversed in the procedure.
There is thus much to be gained from amplifying evanescent waves, provided amplification is limited to spatial frequencies with signal above the noise floor.By equating the signal-to-noise ratio SNR with the amplification factor exp(k z L) for a signal measured at a distance L, we obtain the maximum useful spatial frequency where SNR (in dB) is the signal-to-noise ratio of the amplitude |E x |.A detailed derivation of Eq. ( 4) is presented in the Supplementary Information, accompanied by additional numerical examples in Supplementary Fig. 1.Equation ( 4) directly estimates the maximum spatial frequency that can be retrieved, and thus the resolution that can be achieved, at any z and SNR.Increasing z could thus be advantageous -reducing the near field perturbation induced by the antenna -provided SNR increases accordingly.
1 The transmitted field amplitude is sampled as a function of the time delay of a probe pulse on a photoconductive antenna which probes the radiating near field.The electric field is polarized in x, using the reference frame shown in Fig. 1.Diffraction limited imaging of the finest feature would require a wavelength of λ/2 = d, i.e., frequencies of 0.75 − 1 THz.See Supplementary Fig. 2 for detailed images of samples, near field probe, and experimental layout.
Figure 3(a) and 3(b) shows the measured intensity |E x | 2 emerging from two apertures (diameter and separation: 200 µm) at a frequency of 1.5 THz and 0.38 THz respectively, as well as the associated spatial Fourier transform magnitude | Ẽx |. Figure 3(c) shows the average measured intensity across y = 0 ± 100 µm as a function of x and frequency, and highlights that apertures cannot be discerned directly, either due to diffraction at higher frequencies, or because of evanescent decay at lower frequencies.We calibrate the probe-to-sample distance L by considering the complex field at a frequency above the diffraction limit (here: 1.5 THz), and adjusting L in Eq. (3) to maximize image sharpness (see Supplementary Fig. 3).From | Ẽx |, we then obtain the frequency-dependent SNR via the ratio between the maximum amplitude in the propagating region |k| < k 0 , and the average amplitude in the evanescent region |k| > k 0 , from which we estimate the experimentally accessible k max via Eq.( 4).For the two-aperture experiment we find L = 172 µm (i.e.L 0.87λ at 1.5 THz and L 0.22λ at 0.38 THz) and SNR = 4 − 14 dB between 0.38-1.5 THz, resulting in k max /k 0 = 1 − 1.9 (see Supplementary Fig. 4).We then implement the SL via Eq.( 3), followed by spatial low-pass filtering bounded by k max .Figure 3 Finally, we implement our procedure on a complex, large-area image formed by a laser-written metal sheet containing the letters "THZ" (minimum feature size: d = 150 µm).In this case, the calibration yields L = 440 µm, SNR=15-25 dB, and k max /k 0 = 1 − 3.5 between 0.2-1 THz (see Supplementary Fig. 4). Figure 4 3) at different frequencies with k max = k 0 , i.e., a conventional lens simply reversing phase.At 1.0 THz (i.e., the diffraction limit), the letters "THZ" can be discerned; at 0.67 THz, only the vertical features are resolved; lower frequencies do not provide a sufficiently sharp image to discern the finer features of the original image, with only a single large spot occurring at 0.18 THz. Figure 4(c) shows the corresponding retrieved |E SL x | 2 at different frequencies with k max /k 0 as labelled: vertical features are significantly sharpened.Note that horizontal features do not let E x through: the slits forming the letters act as parallel plate waveguides with width smaller than half a wavelength, in which solely the TEM mode polarized perpendicularly to the thinnest features can propagate, so that thin features in y (x) only appear for the E x (E y ) field.As a result, for x−polarized fields the letters' vertical (y-oriented) features are clearest, in agreement with simulations (see Supplementary Fig. 5).We repeat the above procedure for a polarization oriented in y relative to the sample orientation, and plot the corresponding |E y | 2 in Fig. 4(c) to resolve the horizontal features of each letter.The full image is obtained by summing the two contributions: Fig. 4(d) show the resulting showing the emergence of the letters "THZ" at all frequencies, down to λ/7.These results are in agreement with simulations of the transmitted subwavelength pattern (see Supplementary Fig. 5).Remarkably, the resolution achieved in Fig. 4 is higher than that in Fig. 3, even though the near field antenna's distance to the object was more than doubled, thanks to a higher SNR exceeding the loss from increased evanescent decay.

DISCUSSION
In this paper, we have presented a novel superlensing approach which numerically amplifies measured evanescent fields to obtain complex images with subwavelength features, limited only by the noise of the instrument.We presented experiments illustrating its implementation at THz frequencies using commercially available facilities.Our approach can be adapted to suit any near field experiment which measures amplitude and phase, immediately providing a pathway for increasing the imaging resolution of near field setups at any frequency.High resolution near field measurements in the radiating rather than reactive near field will allow accurate near-field imaging without perturbing the intrinsic field of structures strongly susceptible to local disturbances, such as high-Q/topological resonators [41,42] and photonic crystal defects [43].

Experimental setup
A schematic of the experimental setup is shown in Supplementary Fig. 2. We use a commercially available THz-TDS System (Menlo TERAK15), which relies on THz emission from biased photoconductive antennas that are pumped by fiber-coupled near-infrared pulses (red line; pulse width: 90 fs; wavelength: 1560 nm).Terahertz lenses collimate and focus the beam towards the sample.The THz field emerging from the THzLC is sampled as a function of the time delay of a fiber-coupled probe pulse on another photoconductive antenna, which forms the THz detector.The electric field is polarized in x, using the sample orientation and reference frame shown in Fig. 1.A moveable, fibercoupled near-field (NF) detector module enables the measurement of the x-polarized electric field at the output of the laser-machined samples.The near-field is spatio-temporally resolved at every point via a raster scan (step size: 25-50 µm).Fast Fourier transforms of the temporal response at each pixel position provide the spectral information (spectral resolution: 5-8 GHz).

Noise limit derivation
⊥ , let T be the transfer function of the fields along z: such that Ẽ(k x , k y , z) = T (k x , k y , z) Ẽ(k x , k y , z = 0) (which is propagation in free space) and conversely Ẽ(k x , k y , z = 0) = T (k x , k y , −z) Ẽ(k x , k y , z) (which is the superlensing procedure).For propagating fields k ⊥ < k 0 and For evanescent fields, k ⊥ > k 0 and with a − sign for the evanescent decay (transfer from object to measurement distance) and + sign for the superlensing procedure.
At a distance z from the source, the spatial Fourier transform of the measured field is given by where Ẽm (k x , k y ) is the actual field (i.e., in the absence of noise), and δm (k x , k y ) is the noise due to the measurement instrument, where typically The reconstructed field is given by The first product is exactly the object field we wish to obtain, while the second term is the noise in the reconstructed field δSL = δm exp(|k z z|), where we omit the explicit (k x , k y ) dependence for brevity.Compared to the measurement noise, the (amplitude) noise penalty in dB of the procedure is thus ∆ noise = 10 log 10 δSL δm = 10 log 10 (exp(k z z)) = 10k z z log (10) .
When the noise penalty exceeds the initial signal to noise ratio of the measurement SNR, the image becomes dominated by noise.This occurs when The right hand term increases steadily with k ⊥ , and the limiting value of k ⊥ = k max at which equality is achieved can be obtained by rearranging Eq. (11) as which corresponds to Eq. ( 4) in the main manuscript.Filtering out spatial frequencies higher than k max ensures that the noise penalty due to the superlensing procedure does not generate noise levels above the measurement's signal level, and thus ensures clear images.

Noise limit examples
We now consider an example to showcase the implications of Eq. ( 12) by expanding our analysis on the simulations associated with Fig. 2 of the manuscript.Supplementary Fig. 1(a) (left) shows a simulation of the amplitude E x (x) as a function of normalized propagation length z/λ for the double aperture case of Fig. 2 in the manuscript, before any superlens procedure.Note in particular that inside the interval (2π) −1 < z/λ < 1, i.e., in the radiating near field, the two apertures are not resolved.On the right of Supplementary Fig. 1(a) we show the target amplitude at the source, where both apertures are resolved for all choice of z/λ.The blue line in Supplementary Fig. 1(b) shows the the Fourier transform | Ẽx | at a distance z/L = 0.5, as per Fig. 2(b) of the manuscript.To enable a clear analysis, we add random amplitude and phase, resulting in a nominally flat SNR of 30 dB, shown as a line curve in Supplementary Fig. 1(b).We then perform the superlensing procedure without any filtering, for different propagation lengths z, always adding white noise such that SNR=30 dB before amplification.The resulting | ẼSL x (k x , k y )|, normalized to | ẼSL x (0, 0)|, is shown in Supplementary Fig. 1(c), where the color scale is saturated at unity.For short lengths z/λ, the procedure yields the required high spatial frequency components.For increasingly long propagation lengths however, amplified evanescent wave amplitudes are greater than those of the propagating waves (white regions).Equation 12predicts the boundary between these two regions for SNR=30 dB.To show this, Supplementary Fig. 1(c) also shows contour lines of the function for different choices of ∆ SNR as labelled.For ∆ SNR = SNR = 30 dB (dark blue), where the ∆ SNR value used in Eq. ( 13) is that of the actual SNR, Eq. ( 13) indeed yields the spatial frequency boundary inside which the evanescent field amplitude after the virtual SL remain below that of propagating waves.Decreasing ∆ SNR narrows the available range of k max , as expected, with ∆ SNR = 0 dB corresponding to k max = k 0 .Given a set of experimental conditions, Eq. ( 12) can thus be used as a first estimate of the low-pass spatial frequency boundary to include after the superlensing procedure, and which can be fine-tuned until an suitable image is obtained.This procedure contains a trade-off between image sharpness and noise: A wider k max /k 0 boundary has the effect of deteriorating the image by using amplified noisy spatial frequencies, as shown in Supplementary Fig.Here we use a single feature, highlighted by the white line in (f), and fit a single gaussian.In this case, z = L = 440 µm.In both cases, this numerical processing of our experimental data functionally corresponds to adjusting the focus of a conventional lens, until the sharpest image is obtained.

FIG. 1 .
FIG.1.Concept schematic of virtual superlens.(i) Sub-wavelength spatial features are carried by evanescent waves which exponentially decay over a L (red).(ii) The resulting lower-resolution image is detected by a near field probe.The collected evanescent fields are then numerically amplified over L (green), leading to (iii) the original image, analogously to a superlens (blue).Wavelength-scale information is carried by propagating waves (dashed).
Figures3 and 4showcase our technique on two distinct imaging experiments.Our experiment uses a commercial pulsed THz source (Menlo TERAK15, 0.1-3 THz).Lenses collimate and focus the THz beam towards a patterned laser-machined samples containing sub-wavelength features (d = 150 − 200 µm).The transmitted field amplitude is sampled as a function of the time delay of a probe pulse on a photoconductive antenna which probes the radiating near field.The electric field is polarized in x, using the reference frame shown in Fig.1.Diffraction limited imaging of the finest feature would require a wavelength of λ/2 = d, i.e., frequencies of 0.75 − 1 THz.See Supplementary Fig.2for detailed images of samples, near field probe, and experimental layout.Figure3(a) and 3(b) shows the measured intensity |E x | 2 emerging from two apertures (diameter and separation: 200 µm) at a frequency of 1.5 THz and 0.38 THz respectively, as well as the associated spatial Fourier transform magnitude | Ẽx |.Figure3(c)shows the average measured intensity across y = 0 ± 100 µm as a function of x and frequency, and highlights that apertures cannot be discerned directly, either due to diffraction at higher frequencies, or because of evanescent decay at lower frequencies.We calibrate the probe-to-sample distance L by considering the complex field at a frequency above the diffraction limit (here: 1.5 THz), and adjusting L in Eq. (3) to maximize image sharpness (see Supplementary Fig.3).From | Ẽx |, we then obtain the frequency-dependent SNR via the ratio between the maximum amplitude in the propagating region |k| < k 0 , and the average amplitude in the evanescent region |k| > k 0 , from which we estimate the experimentally accessible k max via Eq.(4).For the two-aperture experiment we find L = 172 µm (i.e.L 0.87λ at 1.5 THz and L 0.22λ at 0.38 THz) and SNR = 4 − 14 dB between 0.38-1.5 THz, resulting in k max /k 0 = 1 − 1.9 (see Supplementary Fig.4).We then implement the SL via Eq.(3), followed by spatial low-pass filtering bounded by k max .Figure3(d) and 3(e) show |E SL x | 2 and | Ẽx | at 1.5 THz and 0.38 THz respectively: the apertures are now resolved.The associated average intensity over the middle of the two apertures (Fig. 3(f)), shows this procedure works over the entire THz frequency band.Finally, we implement our procedure on a complex, large-area image formed by a laser-written metal sheet containing the letters "THZ" (minimum feature size: d = 150 µm).In this case, the calibration yields L = 440 µm, SNR=15-25 dB, and k max /k 0 = 1 − 3.5 between 0.2-1 THz (see Supplementary Fig.4).Figure4(a) shows the measured |E x | 2 at different frequencies as labelled.Figure 4(b) shows the corresponding retrieved field |E SL x | 2 through Eq.(3) at different frequencies with k max = k 0 , i.e., a conventional lens simply reversing phase.At 1.0 THz (i.e., the diffraction limit), the letters "THZ" can be discerned; at 0.67 THz, only the vertical features are resolved; lower frequencies do not provide a sufficiently sharp image to discern the finer features of the original image, with only a single large spot occurring at 0.18 THz.Figure4(c) shows the corresponding retrieved |E SLx | 2 at different frequencies with k max /k 0 Figures3 and 4showcase our technique on two distinct imaging experiments.Our experiment uses a commercial pulsed THz source (Menlo TERAK15, 0.1-3 THz).Lenses collimate and focus the THz beam towards a patterned laser-machined samples containing sub-wavelength features (d = 150 − 200 µm).The transmitted field amplitude is sampled as a function of the time delay of a probe pulse on a photoconductive antenna which probes the radiating near field.The electric field is polarized in x, using the reference frame shown in Fig.1.Diffraction limited imaging of the finest feature would require a wavelength of λ/2 = d, i.e., frequencies of 0.75 − 1 THz.See Supplementary Fig.2for detailed images of samples, near field probe, and experimental layout.Figure3(a) and 3(b) shows the measured intensity |E x | 2 emerging from two apertures (diameter and separation: 200 µm) at a frequency of 1.5 THz and 0.38 THz respectively, as well as the associated spatial Fourier transform magnitude | Ẽx |.Figure3(c)shows the average measured intensity across y = 0 ± 100 µm as a function of x and frequency, and highlights that apertures cannot be discerned directly, either due to diffraction at higher frequencies, or because of evanescent decay at lower frequencies.We calibrate the probe-to-sample distance L by considering the complex field at a frequency above the diffraction limit (here: 1.5 THz), and adjusting L in Eq. (3) to maximize image sharpness (see Supplementary Fig.3).From | Ẽx |, we then obtain the frequency-dependent SNR via the ratio between the maximum amplitude in the propagating region |k| < k 0 , and the average amplitude in the evanescent region |k| > k 0 , from which we estimate the experimentally accessible k max via Eq.(4).For the two-aperture experiment we find L = 172 µm (i.e.L 0.87λ at 1.5 THz and L 0.22λ at 0.38 THz) and SNR = 4 − 14 dB between 0.38-1.5 THz, resulting in k max /k 0 = 1 − 1.9 (see Supplementary Fig.4).We then implement the SL via Eq.(3), followed by spatial low-pass filtering bounded by k max .Figure3(d) and 3(e) show |E SL x | 2 and | Ẽx | at 1.5 THz and 0.38 THz respectively: the apertures are now resolved.The associated average intensity over the middle of the two apertures (Fig. 3(f)), shows this procedure works over the entire THz frequency band.Finally, we implement our procedure on a complex, large-area image formed by a laser-written metal sheet containing the letters "THZ" (minimum feature size: d = 150 µm).In this case, the calibration yields L = 440 µm, SNR=15-25 dB, and k max /k 0 = 1 − 3.5 between 0.2-1 THz (see Supplementary Fig.4).Figure4(a) shows the measured |E x | 2 at different frequencies as labelled.Figure 4(b) shows the corresponding retrieved field |E SL x | 2 through Eq.(3) at different frequencies with k max = k 0 , i.e., a conventional lens simply reversing phase.At 1.0 THz (i.e., the diffraction limit), the letters "THZ" can be discerned; at 0.67 THz, only the vertical features are resolved; lower frequencies do not provide a sufficiently sharp image to discern the finer features of the original image, with only a single large spot occurring at 0.18 THz.Figure4(c) shows the corresponding retrieved |E SLx | 2 at different frequencies with k max /k 0 1(d), left; a narrower k max /k 0 boundary removes high spatial frequencies needed to resolve fine feature sizes, as shown in Supplementary Fig.1(d), right.

Supplementary Fig. 1 . 3 .
Numerical example illustrating the effect of the SNR, when applying the superlens procedure of Eq. (3) in the manuscript.(a) Left: raw simulation of the amplitude Ex(x) as a function of z/λ for the double aperture case of Fig. 2 in the manuscript, before any superlens procedure.Right shows the target amplitude at the source, that is, what a figure should look like after perfect, noiseless superlensing: both apertures are resolved for all choice of z/λ.(b) Dashed blue line shows an example | Ẽx| at a distance z/λ = 0.5, as per Fig. 2(b) of the manuscript.The red line shows the same field after adding random amplitude and phase resulting in a flat SNR of 30 dB.(c) Calculated normalized amplified spatial Fourier transform for the flat SNR = 30 dB, as a function of its normalized propagation length z/λ.Color scale has been saturated to 1 for clarity.Solid curves show Eq.(13) choosing different values of ∆SNR as labelled.Note that the case ∆SNR = SNR = 30 dB corresponds to the boundary between where high spatial frequencies have a comparable magnitude to propagating waves.(d) Resulting E SL x (x) after applying the superlens procedure, using a low-pass filter function bounded by kx/k0 as per Eq. 13 for different values of ∆SNR as labelled.If ∆SNR > SNR, the aperture images are plagued by noise; if ∆SNR ≤ SNR, the maximum distance at which the retrieval procedure produces the image is gradually reduced.The maximum useful z/λ occurs at the point in which kmax/k0 contains the feature needed (here, kmax/k0 = 2.5, which for ∆SNR = 30 dB occurs at z/λ = 0.5.)Retrieval of the sample-to-switch distance z = L. (a) Measured |Ex| at 1.5 THz for the two-aperture experiment, where the spatial frequencies which resolve the aperture features are propagating.(b) Example image when the phase of | Ẽx| is adjusted by −k0z, where z = 200 µm.We consider its field magnitude through the apertures' center (white dotted line), and show it as blue circles in (c).Red curves in (c) show a Gaussian fit to a double Gaussian function, used to obtain the average standard deviation σ.(d) Average fitted σ as a function of z. (e)-(h) Same as (a)-(d), performed for the "THZ" sample.