Terahertz non-label subwavelength imaging with composite photonics-plasmonics structured illumination

Inspired by the capability of structured illumination microscopy in subwavelength imaging, many researchers devoted themselves to investigating this methodology. However, due to the free propagating feature of the traditional structured illumination fields, the resolution can be only improved up to double times compared with the diffractied limited microscopy. Besides, most of the previous studies, relying on incoherent illumination sources, are restricted to fluorescent samples. In this work, a subwavelength nonfluorescent imaging method is proposed based on the terahertz traveling wave and plasmonics illumination. Excited along with a metal grating, the spoof surface plasmons are employed as the plasmonics illumination. When the scattering waves with the SSPs illumination are captured, the high order spatial frequency components of the sample are already encoded into the obtainable low order ones. Then, an algorithm is summarized to shift the modulated SF components to their actual positions in the Fourier domain. In this manner, high order SF components carrying the fine information are introduced to reconstruct the desired imaging, leading to an improvement of the resolution up to 0.12 lambda. Encouragingly, the resolution can be further enhanced by tuning the working frequency of the SSPs. This method holds promise for some important applications in terahertz nonfluorescent microscopy and sample detection with weak scattering.


Introduction
Optical microscopy has become an essential tool in a wide range of applications, including biological imaging, material identification, device creation, etc. However, in conventional far-field microscopy, one of the major limitations is the insufficient resolution, which is constrained to approximately half the working wavelength of the illumination light [1]. This phenomenon is termed the Rayleigh limit or classical diffraction limit [2]. Specifically, the samples' fine information is encompassed in the high-order spatial frequency (SF) components, manifesting as the evanescent waves. They reduce exponentially perpendicular to the propagating direction [3]. Thus, the high-order SF components permanently vanish in the far-field, posing a wavelength level resolution in the ideal circumstance [4].
To overcome this resolution limitation, a series of optical strategies have been developed over the past decades. For instance, the super-resolution fluorescence microscopy techniques, including stimulated emission depletion [5,6], stochastic optical reconstruction [7,8], structured illumination microscopies (SIM) [9][10][11], have endowed a breakthrough in subwavelength resolution imaging. Among the myriad fluorescence super-resolution approaches, SIM is considerably remarkable due to its fast and highly parallelizable feature [12]. In SIM, a finely structured illumination pattern is applied to illuminate the sample. In this case, the high-order SF components (corresponding to the fine detail of the sample) can be down-modulated into the passband of the objective and turn out to be detectable. By post-processing, the highorder SF information can be precisely shifted and located at the true position in the spatial frequency domain (k-domain). As a result, one can defeat the diffraction barrier and finally achieve a subwavelength resolution.
However, owing to the free-propagating feature of the structured illumination fields, the spatial resolution of traditional SIM is restricted to a two-fold enhancement in comparison with the diffraction-limited techniques [13]. Aside from the resolution, the application scenarios of SIM should be taken into consideration. Most of the aforementioned researches depend on incoherent imaging techniques and fluorescent samples, which become out-of-operation for label-free imaging [14][15][16][17][18][19]. From the perspective of the terahertz range, the situation is similar. The existing researches focus on the fluorescent samples and obtain an excellent resolution for subwavelength imaging [20]. However, to the best of our knowledge, the terahertz label-free subwavelength imaging with structured illumination method remains incomplete.
In this work, a terahertz subwavelength imaging method of label-free samples is proposed based on composite photonics-plasmonics structured illumination. The spoof surface plasmons (SSPs) subsisted on a metal grating are employed for the plasmonics illumination. When the sample is impinged by the SSPs, the scattering field is encoded. That is, the high-order SF components of the sample (carrying fine information) are down-transformed to the low-order and obtainable SF ones. Then, the scattering field is caught to get the modulated high-order SF information. Similarly, with vertical and oblique incidences of the traveling waves (photonics illuminations), low-order SF information of the sample is also acquired. Subsequently, an iterative algorithm is applied to extract and locate all the SF components in their true positions. In this way, a wide SF spectrum is generated, which is imperative for subwavelength imaging. Based on this principle, the simulated resolution of the proposed method reaches up to 0.12λ. Note that the resolution can be further improved by elaborately design the SSPs' working frequency. This method is beneficial for applications of terahertz nonfluorescent microscopy and weak scattering sample detection.

Imaging Principle
The imaging principle in this work relies on the coherent scattering from the labelfree samples. Before elaborating on the imaging principle, there are several settings to be clarified. Firstly, this process is between the Fourier domain (k-domain) and the spatial domain (x-y), where k is the SF vector. Here, we consider the condition where the sample is positioned along the x-direction and infinite in the y-direction. The incident wave is on the x-z plane, lacking the component ky. Secondly, we postulate that when an incident wave impinges on a thin subwavelength sample, the scattered wave can be illustrated as s(x) = o(x)exp(iktx). o(x) refers to the sample's complex field [21] and kt is the SF component of the incident wave along the x-axis. Thirdly, the coherent point spread function (PSF) hc(x) of our imaging system is defined via its Fourier transform Hc(kx). Hc(kx) is modeled by an abrupt rectangular low-pass filter, which has a width of 2kc (kc is the cut-off frequency). Since the imaging system of this work is in the free space, kc can be substituted by k0 (the wavenumber of vacuum). Hc(kx) is rendered into gray in Fig.1 (a), where C denotes the module value of the SF component. On this basis, the detected intensity imaging of the sample is indicated as i(x) and derived as the following equation [21] ( ) where  is the convolution operator. Eq. (1) can be written in the k-domain, which is signified with Eq. (2) where I(kx) and S(kx) are the SF distributions of the detected image intensity, and the scattering field, respectively. It should be noticed that is the SF distribution of the sample's field. Eq. (2) thus can be derived as In this manner, the spatial spectrum of the sample is shifted to kx´, manifesting as kx´= kx − kt. The detected image intensity I(kx) contains the attainable SF information, which is governed by |kx´| < kc. Since kc = k0, it is written as −k0 + kt < kx < k0 + kt.
According to the aforementioned principle, the conversion of the collected SF range is attributed to the tangential component kt of the illuminating wave. When a plane wave impinges the sample vertically, the obtained SF region is −k0 < kx < k0 (same as the system's passband), as shown in Fig. 1(a). Analogously, a plane wave with left oblique illumination leads to the acquisition of the SF spectrum in −k0 − k0sinθ < kx < k0 − k0sinθ (θ is the incident angle). For the right oblique incidence, the acquired range is −k0 + k0sinθ < kx < k0 + k0sinθ ( Fig. 1(b)). Since θ < 90°, the detectable region is limited by ±2k0.
To further increase the observable regime, the evanescent wave is exploited as the illuminated wave (along both x and −x-directions). The tangential component keva of the evanescent wave is larger than k0. Consequently, the attainable SF regime (−k0 ± keva < kx < k0 ± keva) is extended further, as depicted in Fig. 1(c).
Finally, through judiciously tailoring the incident angle of the plane wave and the tangential component of the evanescent wave, a wide integrated k-domain information can be spliced with no empty regions. It is noted that the wide information in the kdomain is imperative for high-resolution imaging. Thus, applying the inverse Fourier transform, the subwavelength spatial image of the sample is restored.

Reconstrcution algorithm
After gathering all of the raw intensity images under different illuminated patterns, a meliorated Gerchberg-Saxton algorithm [22,23] is exploited to iteratively generate a wide k spectrum. In this way, a subwavelength-resolution image can be realized. The procedure involves five steps.
Step 1. An initial value Is 1/2 e iφ of the subwavelength sample in the spatial domain is assumed, using φ = 0 and Is being a constant value (detected data group is also available [1]). Applying the Fourier transform for the initial value, a wide band spectrum W(kx) can be presented in the k-domain Step 2. The low SF components are restored with the image under vertical illumination, i.e., oblique illumination with an incident angle of 0°. A subregion of the wide k-domain is extracted and dealt with the invert Fourier transform, producing a new temporary image It 1/2 e iφₜ The selected area in the k-domain hinges on the low-pass filter Hc(kx). The location of the low-pass filter depends on the tangential wavenumber of the illumination pattern. Under the vertical incidence (tangential wavenumber being zeros), the position of the filter's center is kx = 0. Thus, the domain of this subregion is −k0< kx < k0. Resultantly, the aforementioned procedure is explicated as where F −1 indicates the inverse Fourier transform. Then, the amplitude It 1/2 of the temporary image It 1/2 e iφₜ is substituted by the detected image amplitude Id 1/2 . Next, we use the Fourier transform of Id 1/2 e iφₜ to replace the corresponding subregion in the kdomain. An updated SF spectrum Wnew(kx) is configured Step 3. By repeating step 2 (select a small, circular region of k-domain and update it by the corresponding detected data), the high SF components are also obtained with oblique and evanescent illuminations. For the left and right oblique illumination with incidence angle θ, the extracted subregion of k-domain is displayed as −k0 − k0sinθ < kx < k0 − k0sinθ and −k0 + k0sinθ < kx < k0 + k0sinθ, respectively. Similarly, the extracted subdomain with the evanescent illumination is presented as −k0 ± keva < kx < k0 ± keva, where keva is the evanescent wave's component along the x-direction.
Step 4. Repeating steps 2 and 3 until a self-consistent solution is realized.
Step 5. The converged solution Wsolved(kx) is processed with an inverse Fourier transformation to recover a subwavelength image isolved(x). After an origin calibration, a subwavelength-resolution image is produced. The complete pseudocode is demonstrated in the supplementary material.

SSPs Illumination Field
As mentioned above, evanescent waves are applied as the illumination field to achieve a down-modulation of the high SF information. In the optical region, the surface plasmon polaritons (SPPs), produced by collective oscillations of electron gases near the metallic surfaces and interfaces [24], can be utilized as evanescent wave illumination [1,[25][26]. However, metals display the ideal electric conductor property in the terahertz zone, as a result of which natural SPPs are severely declined [27]. To solve this problem, the evanescent wave illumination is achieved by the SSPs in this work. The SSPs propagate on a metal grating surface and decay exponentially along the vertical direction [28,29]. Via manipulating the working frequencies of the SSPs, a wealth range of tangential wavenumbers (keva) are available, benefiting the system design.
The schematic diagram of the proposed approach is shown in Fig. 2, wherein a metal grating is applied to sustain the SSPs on the surface. Two thin samples with different relative permittivities are deposited above the grating in a subwavelength distance. Leveraging on a coupler between the SSPs and traveling waves, the SSPs are excited adequately in this device. Then, the SSPs propagate along the grating and impinge on the samples, rendering a detectably scattering wave. After that, the SSPs propagate forward and exit as a traveling wave through another coupler. It is worth emphasizing that the SSPs can transmit along both the x-direction and −x-direction, just depending on which coupler it is excited on.
Specifically, the grating's unit cell is presented in Fig. 3(a), where h = 0.049 mm, p = 0.02 mm, and d = 0.01 mm (the same size in all simulations). The dispersion of the SSPs contains the fundamental mode and the high-order modes [30]. To prevent interference from multiple modes, only the fundamental mode [ Fig. 3(a)] is employed in this work. In this work, the working frequencies of SSPs are f1 = 1.1 THz, f2 = 1.23 THz and f3 = 1.33 THz. Corresponding to the working frequencies, the tangential wavenumbers (x-direction) are keva1 = 1.387k0, keva2 = 2.082k0, keva3 = 3.539k0. k0 is caculated as 2f0/c, where c is the velocity of light. f0 is the working frequency of the . Notably, further promotion of the imaging resolution is feasible by adding new SSPs illumination with high working frequency. With a higher working frequency, larger SF components can be down-modulated into the passband of the system. Thus, it leads to a larger effective numerical aperture NA = (max{keva1, keva2, keva3, keva4, … } + max{feva1, feva2, feva3, feva4, … } / f0)k0 and a higher resolution [calculated as λ0/(2NA)]. Since more illumination patterns causes a more sophisticated post-process, a balance between the resolution and number of illumination patterns should be elaborately determined.

Simulation and Results
Hereto, the procedure of how to recover a subwavelength image in the terahertz region is perspicuous. To encode the Fourier spectrum, there are various illumination patterns, including vertical, oblique, and evanescent illuminations. For the data collection, the scattering field intensity of the sample is directly detected, which avoids extra laser to excite the sample [20]. For the image reconstruction, a post-process algorithm is applied.
With the assistance of the commercial software COMSOL Multiphysics 5.4a, the scattering fields intensity are investigated under diversiform illuminations. Take the single sample (0.15λ0×0.1λ0 on the x-z plane) imaging as an example, for the vertical illumination, the sample is directly impinged by a plane wave [see Fig. 4(a)]. Similarly, plotted by Figs. 4(b) and (c), the oblique patterns employ plane waves with ±66° incident angles. The corresponding tangential wavenumber are ±0.914 k0. For the evanescent illumination, the sample is placed above the grating surface with an extremely close distance, which is approximately 0.024λ0−0.12λ0. Figs. 4(d) and (e) demonstrate the scattering fields intensity with SSPs propagating along x and -xdirections (working frequency is f1). In addition, the SSPs illuminations at f2 and f3 are also simulated. Consequently, the tangential wavenumbers of the nine illumination waves are [0, 0.914k0, −0.914k0, 1.387k0, −1.387k0, 2.082k0, −2.082k0, 3.539k0,  Fig. 4(f), the nine obtained SF spectrums are presented. Since the working frequencies of propagating waves and SSPs are distinct, the corresponding half-width of the detectable SF regions are derived as (fi / f0) × k0 = 2πc / fi, i = 0, 1, 2, and 3.
Next, the scattering fields intensity under these nine illumination patterns are collected, without the procedure of laser excitation. For the resolution investigation, the length of the detective line on the x-z plane is 2540 μm, which is 268 μm away from the sample. With the post-process algorithm, a wide k spectrum is reconstructed. Then the k spectrum is processed by the apodization and zero-padding function to extract the critical SF information and reveal smooth results, respectively.
For the recovered results, the first demonstration is the imaging of a single, isolated sample with the subwavelength feature, whose dimension is 0.15λ0 × 0.1λ0 on the x-z section. The relative permittivity of the sample is ε = 2.5, characterizing the dielectric property. The reconstructed images of this sample are depicted in Fig. 5(a), which contain various illumination groups. The black curve illustrates the result generated by the whole nine illuminations (group1), presenting a 0.11λ0 full width at half maximum (FWHM). When the SSPs illuminations at f2 and f3 are eliminated in the reconstruction process (group2), the FWHM becomes larger than before, as designated by the blue curve. Ulteriorly, if only the vertical and oblique illuminations (group3) remain without any SSPs illumination, the FWHM increases further. The modulus of the k spectrum is plotted in Fig. 5(b), corresponding to these assorted illumination combinations. The orange region represents the attainable SF spectrum with vertical and oblique patterns. The blue plus orange regions are the whole achieved SF information under the group2 illumination. When the combination is group1, the homologous SF region includes the blue, orange, and gray pieces. Apparently, high SF information in the k-domain is removed without the relevant SSPs illumination, leading to a wide FWHM. From the point of Abbe's criterion (imaging a single sample) [31], this phenomenon is rational. According to the criterion, smaller FWHM means higher resolution, which is attained through collecting more fine information (high SF regions). Since the SSPs illumination brings high SF in k-domain, abandoning them naturally causes the low resolution and large FWHM. Finally, based on Abbe's criterion, this work can realize a 0.11λ0 resolution, which can be further improved via adding SSPs illumination of large keva.
Besides, to verify the robustness of this method, the sample is repositioned along the x-axis. Assuming the sample is placed at the origin for the aforementioned discussion. Whereupon we move it to ±0.5λ0 positions, keeping the other conditions fixed. The result is displayed in Fig. 5(c), which testifies the method's stability with disparate locations. There exists a slight difference between the two curves, which is mainly caused by the personalized iteration number when a convergent solution is generated. Another issue needed to be explored is whether the sample's dielectric feature has an impact on the reconstruction. The relative permittivity is set as 1.5, 2.5, 3.5, and 4.5 with the other conditions being identical. Comfortingly, Fig. 5(d) illustrates that the restored imaging keeps stable, meaning that the method is robust to the variation of the relative permittivity within a certain range.
Up to here, we have demonstrated the 0.11λ0 resolution and the robustness of this method. However, there are many other criteria to evaluate the imaging resolution in the literature. The most likely accepted is the Rayleigh criterion [2], which confirms the minimum resolvable distance by the imaging of two closely separated samples. In the following, the examination of the resolution in the view of the Rayleigh criterion is presented.
First, we consider two same samples with ε = 2.5, h = 0.1λ0, and width ws. They are separated by a gap, whose length is wg = ws. We recover the images with ws ranging from 0.15λ0 to 0.225λ0. Figs. 6(a)-(c) depict the images under these conditions, which confirm the method's subwavelength imaging ability of two samples. Furthermore, in Fig. 6(d), the samples are set extremely small with h = ws = 0.02λ0 for the judgment of the highest resolution. The distance between the center of them is wc = 0.12λ0. In this case, the samples can be approximately regarded as two points with a separate distance approximately being 0.12λ0. To increase the interaction between the tiny samples and the illuminating waves, their relative permittivities are raised to 12. Under this condition, the result verifies that a resolution of 0.12λ0 can be achieved. Subsequently, when the samples get closer to each other (0.11λ0), the image becomes inexplicit. As plotted in Fig. 6(d), one peak is almost a quarter lower than the other, illustrating the method cannot distinguish two samples separated by 0.11λ0. Thus, on basis of the Rayleigh criterion, the ultimate resolution of this method is 0.12λ0. Actually, there exists a slight deviation about the resolution with the Abbe and Rayleigh criterion. One reason is that the samples in the Rayleigh criterion has a 0.02λ0 width, and the gap between them is actually shorter than 0.12λ0. Besides, the undesired reflection between the two samples convey noise to the reconstruction and lead to a slight deviation.
Hereto, the theoretical resolution and the simulated resolution can be exhibited and compared. For this work, the effective NA is decided by the maximum keva of the SSPs illumination, which is 3.539k0. Thus, the highest detectable SF in k-domain is ±(3.539 + f3/ f0) k0 = ±4.66k0, i.e., NA = 4.66k0. In this way, the theoretical resolution is approximately 0.107λ0. According to Abbe's criterion, the simualted resolution is 0.11λ0. Considering the approximation error, it is almost consistent with the theoretical value. With the Rayleigh criterion, the established resolution is 0.12λ0, which has a slight excursion. The offset may caused by the finite size of the samples and the undesired reflection between them. Finally, the imaging resolution is demonstrated as 0.12λ0, almost four times that of diffraction-limited imaging.
Besides, we also restore images of two samples with various relative permittivities. At the beginning, we employ two samples (h = 0.1λ0, ws = 0.15λ0, wg = 0.2λ0) with the same relative permittivity ε = 2.5 as a reference standard, presented in Fig. 7(a). Afterward, the relative permittivities of the left and right samples are stated as ε1 and ε2, respectively. In Fig. 7(b), it is displayed that when ε1 = 2.5 and ε2 = 4.5, the right peak's height is larger than the left one. Analogously, for ε1 = 2.5 and ε2 = 1.5, the left peak is higher than the right peak, shown in Fig. 7(c). Actually, due to the sample's different scattering power with diverse relative permittivities, there emerges unequal peaks. Intuitively, a large relative permittivity yields strong scattering under the impinging wave. To some extent, this method can distinguish two samples with distinct relative permittivities, which may be helpful for the subwavelength dielectric determination. It should be mentioned that the ambiguous image in Fig. 6(d) also has unequal peaks, demonstrating the resolution limit. This means when one peak behaves lower than the other, the reason is not only the inconsistent ε. Thus, this function of relative permittivity distinction can be exploited to two same-size samples with relatively large interval.
Finally, the imaging ability of three samples is also demonstrated in Fig. 7(d). Under this condition, the detected line is moderately adjusted and extended to fit the extended space of the samples. There are three same samples (h = 0.1λ0, ws = 0.15λ0, ε = 2.5) with the gap distances ws1 = 0.25λ0 and ws2 = 0.45λ0. Consequently, they are successfully reconstructed in the recovery image.

Discussion
Firstly, the resolution can be further improved via adding new SSPs illumination of high frequency. In this manner, the effective NA is increased so that the resolution can be enhanced. Remarkably, the working frequencies need to be advisably chosen to ensure a wide k spectrum without interrupted region [22]. Secondly, the sample's scattering field intensity is collected straightly, circumventing an extra laser excitation process. Thirdly, the noise in this system may cause side lobes and artifacts in the reconstructed imaging, which can be further reduced through an optimization of the algorithm [32]. Besides, the influence of couplers' reflection and the higher harmonics on the grating can also be sufficiently suppressed to improve the imaging performance. Fourthly, the method can be exploited to the two-dimensional imaging, the key step of which is to realize multidirectional SSPs illuminations on the x-y plane. One can rotate the grating around the z-axis and simultaneously fix the sample to acquire the multidirectional SSPs. In addition, the multidirectional SSPs illuminations may also be achievable based on an integrated multi-port waveguide [33,34].

Conclusion
In conclusion, this work presents a theoretical proposal to realize a subwavelength coherent imaging method for non-label samples in the terahertz region. A tunable composite photonics-plasmonics illumination pattern is used in the scheme. Depending on the SSPs on a metal grating, the plasmonics illumination pattern is obtained. When the SSPs impinge on the sample,the high-order SF information is down-modulated and contained in the detected scattering field intensity. Accordingly, under the vertical and oblique propagating waves (photonics) illumination, the relatively low-order SF information is also captured. In this manner, an integral wide Fourier spectrum of the sample is achieved, then a subwavelength recovery image is achieved by a post-process algorithm. With nine illumination patterns, a resolution up to 0.12λ0 is verified against both Abbe's and Rayleigh's criteria. Besides, the robustness of this proof-of-concept method for different samples' parameters is analyzed. Finally, improvements on the resolution and image quality are also discussed for further work. We believe that the proposed method offers the potential foundation for terahertz non-fluorescent microscopy and sample detection with weak scattering.