Object-dependent spatial resolution of the reflection-mode terahertz solid immersion microscopy

Terahertz (THz) solid immersion microscopy is a novel promising THz imaging modality that overcomes the Abbe diffraction limit. In our prior work, an original reflection-mode THz solid immersion microscope system with the resolution of 0.15λ (in free space) was demonstrated and used for imaging of soft biological tissues. In this paper, a numerical analysis, using the finite-difference time-domain technique, and an experimental study, using a set of objects with distinct refractive indexes, were performed in order to uncover, for the first time, the object-dependent spatial resolution of the THz solid immersion microscopy. Our findings revealed that the system resolution remains strongly sub-wavelength 0.15–0.4λ for the wide range of sample refractive indices n= 1.0–5.0 and absorption coefficients α= 0–400 cm−1 (by power). Considering these findings, two distinct regimes of the THz solid immersion microscopy were identified. First is the total internal reflection regime that takes place when the sample refractive index is relatively low, while the sub-wavelength resolution is enabled by both the evanescent and ordinary reflected waves at the interface between a high-refractive-index material and an imaged object. Second is the ordinary reflection regime that occurs when the sample refractive index is high enough, so that there is no more total internal reflection at the interface, while only the ordinary reflected waves inside a high-refractive-index material are responsible for the sub-wavelength resolution. The resultant conclusions are general and can be applied for analysis of solid immersion lenses operating in other spectral ranges, such as visible and infrared, given linear nature of the Maxwell’s equations. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Solid immersion lens (SIL) was introduced in 1990 to achieve the sub-wavelength-resolution in the visible microscopy [1]. Since then, it was vigorously explored, transferred to other spectral ranges, and found a variety of applications in different branches of science and technology. Solid immersion effect is a reduction in the electromagnetic beam caustic (focal spot) dimensions, when the beam is focused in free space at a small distance from a high-refractive-index SIL lens, and is due to contribution of evanescent waves from the total internal reflection (TIR) effect [2]. The SIL near-field imaging improves significantly spatial resolution (by a factor of the lens refractive index) and, in certain cases, overcomes the Abbe diffraction limit. This imaging modality retains high energy efficiency due to the lack of sub-wavelength apertures or probes in its optical scheme.
So far, visible and infrared SIL-based optical systems were applied in optical data storage [3][4][5], thermal [6,7] and Raman [8] imaging, photolitography [9], semiconductor science and nanophotonics [10][11][12][13], as well as non-destructive evaluations of electronic circuits [14,15]. Most recently, the SIL imaging was transferred to the terahertz (THz) range [16,17], where our original THz solid immersion microscope arrangement demonstrated the spatial resolution of 0.15λ (λ is a wavelength in free space) and was adapted for reflection-mode imaging of soft biological tissues [18]. The THz SIL technique was recently applied for imaging of printed circuit boards and other dielectric objects [19], as well as soft biological tissues of different nature, such as blades of leaves, tissue spheroids for bioprinting applications, and tissues of the human breast and tongue ex vivo [18][19][20][21][22]; all these measurements revealed strongly inhomogeneous character of tissues at the scale posed by the THz wavelengths.
Several papers were dedicated to theoretical and numerical analysis of electromagnetic beam focusing using SIL [19,[23][24][25][26][27]. Different analytical and computational methods of electrodynamics were applied to study different SIL arrangements, electromagnetic beam polarization, apodization of beam aperture and wavefront aberrations, interplay between evanescent and transmitted waves, and even manufacturing errors and misalignments of SIL components. Despite such a comprehensive theoretical background, a problem of the object-dependent SIL resolution, which was identified in a number of earlier works [17,18,24], and which is of considerable importance for practical applications, still has not been addressed properly. More generally, performance of any near-field optical system (including SIL), regardless of the spectral range and the operation principle, should depend on the optical properties of an imaged object.
For the SIL microscopy, the object-dependent resolution ]is due to variable TIR conditions at the interface between the high-refractive-index SIL lens and a sample. This also results in the object-dependent beam spot formed at the shadow side of a SIL via interference of the propagating and TIR (evanescent in the object) waves. While other modalities of the near-field scanning probe optical microscopy are available, they rely on different physical principles than SIL [28,29]. Many of such methods provide a sub-wavelength spatial resolution thanks to a strong confinement of electromagnetic waves at the scanning probe, such as a sub-wavelength-scale aperture, a dielectric fiber with a strong confinement of guided modes, or a metal cantilever/tip [30][31][32][33][34][35]. For all these modalities, one can also expect the object-dependent resolution, since the conditions of electromagnetic-wave scattering on a sub-wavelength probe should be sensitive to the optical properties and geometrical parameters of a probed object. Evidently, analysis of the object-dependent performance of any near-field imaging technique should be performed by considering the underlying physical principles and peculiarities of an imaging system; while in this study we focus specifically on the object-dependent performance of the reflection-mode SIL microscopy.
The main goal of this paper is to quantify the impact of the object optical properties on the SIL resolution. To this end, we conduct numerical and experimental studies of the object-dependent spatial resolution using an original THz SIL arrangement, reported in Refs. [18][19][20], that features 0.15λ resolution in free space. Figure 1 shows schematic of our THz SIL lens, which is comprised of three main components: • a rigidly-fixed wide aperture aspherical singlet (lens) [36,37] with the focal length of 15 mm and the entrance pupil diameter of 25 mm, made of High-Density Polyethylene (HDPE); • a rigidly-fixed hypohemispherical lens with the diameter of D = 10 mm and the thickness of l = 4.65 mm, made of High-Resistivity Float-Zone silicon (HRFZ-Si); the lens is mounted in front of the focal plane in such way that its spherical surface is concentric with the convergent THz wavefront, while its planar surface is perpendicular to the optical axis [18,19]; • a movable HRFZ-Si window with the thickness of l ′ = 0.25 mm and the diameter of ≃ 50 mm; the window is in contact with the flat surface of the HRFZ-Si hemisphere [18,19]. The distance between the HDPE aspherical lens and the HRFZ-Si hemisphere is z 0 = 2.7 mm. The HRFZ-Si window is mounted on a motorized 2D linear translation stage, which allows displacement in the lateral plane and raster scanning of an imaged object with the THz beam caustic [18]. In the THz range, HRFZ-Si features the high refractive index of n Si ≃ 3.415, as well as negligible absorption and chromatic dispersion [38]. The maximal aperture angle of our THz SIL is θ max = 40 • .
In our numerical analysis, we use the Finite-Difference Time-Domain (FDTD) technique for studying beam propagation in a SIL system [39], while experiments involve imaging of a set of objects with distinct optical properties. The observed numerical and experimental data agree well and demonstrate that the analyzed THz SIL provides strongly sub-wavelength resolution 0.15-0.4λ, within the considered ranges of sample refractive indices n = 1.0-5.0 and absorption coefficients α = 0-400 cm −1 (by power). As a results, we distinguish two regimes of the SIL operation. First is the TIR regime that occurs when the sample refractive index is relatively low, while the sub-wavelength resolution is caused by both the evanescent and ordinary reflected waves at the interface between a high-refractive-index material and a sample. Second is the ordinary reflection regime that occurs when the sample refractive index is high enough, so that there is no more TIR effect and, thus, evanescent waves at the interface are not excited, while only the ordinary reflected wave inside a high refractive index material is responsible for the sub-wavelength resolution. The resultant conclusions have general character and can be applied to analyze performance of SIL systems from other spectral range given linear nature of the Maxwell's equations [40].

Origin of the object-dependent THz SIL resolution
We now consider the origins of high spatial resolution of the reflection-mode SIL-based imaging systems, where an interplay between the evanescent and ordinary reflected waves is observed at the shadow side of a high-refractive-index material [2,18], as illustrated in Fig. 2(a).
Evanescent waves propagate along the SIL flat surface, while their amplitude decays exponentially with the distance from SIL/sample interface [41]. The radiation energy transfer occurs Numerical analysis of the THz SIL resolution at λ = 500 µm when imaging objects with different refractive indexes n obj and absorption coefficients α obj (by power): (a) schematic of the THz SIL with an object placed behind the HRFZ-Si window; (b) field intensity I (r) near the focal point for the free space operation (n obj = 1.0, α obj = 0); (c) normalized resolution δ = FWHM OY /λ as a function of n obj , α obj ; (d)-(f) focal spot shapes for the representative values of the object optical properties: n obj = 1.2, α obj = 0 correspond to the TIR regime of the SIL operation; n obj = 2.4, α obj = 0 and n obj = 3.7, α obj = 0 correspond to the transmitted-wave regime with negative (n Si >n obj ) and positive (n Si <n obj ) dielectric contrast, respectively. In (c), red markers define a set of points, measured experimentally in section 4 (see Table 1, Figs. 3 and 4); while green dashed line shows a cross-section of the color map to be compared with the experimental data (see Fig. 4).
solely along the interface, and, thus, the evanescent field wave vector is directed along the interface k evan = n Si k 0 sin θ i , where k 0 is a wave vector in free space, and θ i is an incidence angle (see Fig. 2(a)). Therefore, k evan is a function of the HRFZ-Si refractive index n Si , which also reduces the effective wavelength of the evanescent field by a factor of n Si sin (θ i ). In turn, for the ordinary reflected wave, which propagates inside HRFZ-Si lens and reflects from the HRFZ-Si / object interface, the wavelength is also reduced by the factor of n Si similarly to conventional liquid immersion microscopy. Evidently, resolution associated with both the evanescent and ordinary reflected waves is enhanced thanks to the high value of the SIL refractive index n Si .
Next, we notice that the ratio between powers carried by the evanescent and ordinary reflected waves of the THz beam changes when varying the object refractive index n obj , as the value of the critical angle θ TIR is given by sin θ TIR = n obj /n Si (see Fig. 2(b)). Particularly, when the object refractive index is high, so that the TIR angle θ TIR become larger than the beam aperture θ max , the TIR part of the THz beam disappears, while only the ordinary reflected waves are responsible for the system resolution. Such operation regime is quite similar to the ordinary immersion microscopy. At the same time, when refractive index of an object is small, so that the TIR effect occurs at the interface, both the evanescent and ordinary reflected waves are responsive for the high image resolution. However, accurate analytical estimation of a resultant THz focal spot geometry, which is formed by these two interfering waves seems to be a daunting task. Therefore, for a more quantitative analysis, we resort to numerical modeling.

Numerical analysis of the object-dependent THz SIL resolution
Dependence of the THz SIL resolution on the sample optical properties was studied using the FDTD method [39]. In Fig. 2(a), a plane wave with the free space wavelength of λ = 500 µm (the frequency of ν ≃ 0.6 THz) was used to radiate optical system from the left side. The classic Yee's approach was used for the spatial-and time-domain discretization of the simulated volume [42]. The Total-Field / Scattered-Field (TFSF) method was used for introducing an ideal plane electromagnetic wave into the simulated volume [43]. The 2 nd -order Mur's absorption boundary conditions allowed to prevent unphysical back-scattering of waves from boundaries of the simulated volume [44]. The spatial step was chosen to be much smaller that the electromagnetic wavelengths (both in free space and in the materials) ∆ r ≃ λ/100, while the time-domain step satisfies the Courant-Friedrichs-Lewy condition ∆ t ≤ ∆ r /(c 0 √ 3), where c 0 = 3 × 10 8 m/s is the speed of light in free space [39].
Although details of the optical system performance differ in the 3D and 2D (cylindrical lens) cases [45], here, we use 2D geometry over 3D in order to speed up simulations and explore wider parameter range, while still capturing the essence of the analyzed effects. The results of 2D simulations are expected to be similar to those of a 3D case up to a multiplier, thus, providing a good first approximation for the resolution analysis of our THz SIL. Particularly, the optical system under study is assumed to be comprised of cylindrical components with the transverse profiles (see Fig. 2(a)) equal to those of a 3D case (see Fig. 1). In 2D FDTD simulations, polarization of the scattered THz field can be either Transverse Magnetic (TM) with the electric field vector E perpendicular to the simulation plane or the Transverse Electric (TE) with E in the simulation plane [39]. In prior numerical studies of our THz SIL system, we found that TM and TE polarization states give very similar results, while the experimentally measured resolution of our system (in free space) was almost polarization independent [17][18][19]; therefore, here, we consider only the TM polarization state for simplicity.
In Fig. 2(a), we show schematic of the modeled THz SIL, behind which an imaged object with refractive index n obj and absorbtion coefficient α obj (by power) is placed. We model the object as a semi-infinite homogeneous rectangular medium in order to eliminate impact of the back-scattered THz waves (from the fascets of the finite-size object) on the THz beam caustic formation at the focal point.
In Fig. 2(b), a typical example of the electromagnetic field intensity distribution I (r) is shown in case of a free space behind the HRFZ-Si window (n obj = 1.0, α obj = 0). It is calculated as where r is a radius vector, t 0 and T is the starting time and the integration interval, respectively, and E (r, t) is an electric field vector. In Eq. (1), the integration is started after electromagnetic wave propagates and fills the entire simulation volume, which leads to the stationary distribution I (r); T equals to the duration of a few full cycles of electric field oscillations T = mλ/c, where m ∈ N.
In Fig. 2(c), relative resolution δ of THz SIL is presented as a function of the sample optical properties n obj , α obj , where δ is estimated as the Full-Width at Half Maximum (FWHM) of the focal spot and normalized by the wavelength λ = 500 µm. For the considered THz SIL with the maximal aperture (half angle) θ max = 40 • , TIR effect at the HRFZ-Si -object interface occurs when the refractive index of an object is n obj = n Si sin (θ max ) ≤ 2.2. Therefore, in Fig. 2(c), two distinct modes of our THz SIL operation are identified: i) the TIR regime (1.0 ≤ n obj ≤ 2.2), with both the evanescent and ordinary reflected waves excited at the HRFZ-Si / object interface; ii) the ordinary reflection regime (n obj >2.2), for which no evanescent modes are excited at the interface and only the ordinary reflected waves exist.
Moreover, within the ordinary reflection regime, two more modalities can be distinguished: ii.a) a negative dielectric contrast (2.2 ≤ n obj ≤ n Si ) at the HRFZ-Si / object interface; ii.b) a positive dielectric contrast (n Si <n obj ).
In Figs. 2(d)-(f), normalized focal spots are shown for each of the abovementioned regimes. From Fig. 2(c) we notice that, although the spatial resolution δ of the THz SIL generally decreases with increasing refractive index n obj and power absorption coefficient α obj of an imaged object in all operation regimes, it remains strongly sub-wavelength withing the considered ranges of the sample optical properties variation 1.0 ≤ n obj ≤ 5.0 and 0 ≤ α obj ≤ 400 cm −1 . We also notice an oscillating character of the spatial resolution as a function of n obj at α obj ≤ 150-200 cm −1 , which can be attributed to the THz wave interference inside the HRFZ-Si hemisphere. The highest resolution of δ ≃ 0.15 is observed for the objects with relatively low refractive indices n obj ≤ 2.7 and absorption coefficients α ≤ 50 cm −1 . Nevertheless, our numerical analysis predicts high resolution δ ≤ 0.2-0.25 even when imaging highly absorbing media, or at local minima in Fig. 2(c).

Experimental study of the object-dependent THz SIL resolution
In our experimental studies, we used a set of test objects fabricated from different optical materials having distinct refractive indexes n obj and absorption coefficients α obj (by power); among them: Polymethylpentene (TPX), Polytetrafluoroethylene (PTFE or Teflon), High-Density PolyEthylene (HDPE), magnesium fluoride (MgF 2 ), crystalline SiO 2 , zinc selenide (ZnSe), gallium phosphide (GaP), sapphire (α-Al 2 O 3 ), and germanium (Ge). In Tab. 1 and Fig. 2(c), optical properties of these materials at λ = 500 µm are summarized, with refractive indices for both ordinary n o and extraordinary n e rays identified.
The test objects with an abrupt (step-like) change in the reflectivity, were fabricated as follows. Rectangular blocks of the test materials were places in contact with polished metal (Al) plates and, then, submerged into epoxy resin to form a monolithic compound. Together with epoxy, the Table 1. THz optical properties n obj , α obj of the experimentally studied objects at the free-space wavelength λ = 500µm. For anisotropic crystalline media, refractive index is presented for both ordinary n o and extraordinary n e rays. sample were then grinded and polished in order to obtain a flat surface with an abrupt metal / test material interfaces, as shown in Fig. 3(a). For the THz imaging, a continuous-wave THz SIL microscope operating in the reflection mode and detailed in Refs. [18,19] was used. THz images are recorded by raster scanning the test objects with a focused THz beam. Backward-wave oscillator [53] is used as a THz-wave emitter at λ = 500 µm (ν ≃ 0.6 THz), and a Golay cell (an opto-acoustic cell) [54] is used as a THz power detector. Furthermore, a 22-Hz-mechanical chopper is used to modulate the THz beam, as required for the Golay cell operation. Moreover, lock-in detection is used to improve the signal-to-noise ratio in the THz images. An attenuator is used to prevent the detector overload. A Keplerian telescopic system with a 1-mm-diameter metal diaphragm in the focal plane is used to homogenize the THz-beam intensity over its aperture.
In Fig. 3(b)-(k), results of the test-object imaging at λ = 500 µm using our THz SIL microscope are shown. For each of the test objects, we show a step-like reflectivity profile, related crosssection of the THz image intensity I im (y) = I im (x, y) | x=x ′ (here, x ′ is some fixed coordinate), as well as estimation of the THz focal spot profile I fs (y) are presented. The THz focal spot profile I fs is calculated as a modulus of the first derivative of the intensity profile I im [18,19] I fs (y) ∝ Finally, resolution of the analyzed THz SIL microscope was estimated experimentally as a FWHM of the first maximum in thus calculated THz focal spot profiles. One can notice asymmetric character of the THz focal spots, as well as strong side lobes in their intensity distributions observed experimentally. This is a typical effect in SIL systems, which appears to be due to the asymmetric excitation of the evanescent an reflected waves at the high-refractive-index material / object interface driven by inhomogeneous structure of an imaging object [17].
In Fig. 4, we present experimentally estimated spatial resolution, which is shown as a function of a sample refractive index n obj , as well as vertical and horizontal error bars. The vertical error bars account for the one sigma confidence interval of the resolution estimation. Particularly, for each test object, we calculated the average resolution and an associated error by averaging over several scans across a material/metal interface performed at different locations along the interface. In turn, the horizontal error bars define possible uncertainties in the value of the refractive index for anisotropic materials (see Tab. 1) due to generally unknown crystallographic orientation of such materials with respect to the incident plane wave of a source. The experimental results are overlapped with the corresponding 2D-FDTD numerical data from Fig. 2(c) and with the results of a supplementary numerical analysis of our THz SIL system, carried out using a quasi-3D Finite-Element Frequency-Domain (FEFD) method [55] within the COMSOL Multiphysics Software and aimed at uncovering the differences between 2D and 3D representations of a SIL systems.

Fig. 4.
A comparison between the experimentally estimated (red circles) and numerically predicted (using 2D-FDTD and quasi-3D-FEFD; see green squares and orange area, respectively) resolution parameter δ of the THz SIL microscope for the considered set of the test objects (see Table 1), with their refractive indices n obj and absorbtion coefficients α obj at λ = 500 µm. Results of the 2D-FDTD simulations are taken from Fig. 2(c). Data of the quasi-3D-FEFD analysis consider variations of the THz-beam spot geometry and, thus, the resolution parameter δ with respect to the orientation of the electric field vector E in the object plane; while the resolution is better in the direction perpendicular to the electric field E. Vertical error bars of the experimental data correspond to the one sigma confidence interval of resolution estimation, while horizontal ones correspond to possible uncertainties in the refractive index while measuring anisotropic materials. Experimental markers are numbered according to the list of test materials fromin Table 1, while experimental resolution for the THz SIL operation in free space (n obj = 1.0) is from Refs. [18,19].
A quasi-3D FEFD method models optical systems featuring a rotational symmetry with respect to the set optical axis (OZ in our case). It is carried out using an optical system 2D cross-section in cylindrical coordinates (r, z); where r is a radius vector perpendicular to OZ axis. Cylindrical Perfectly-Matched-Layers (PML) were applied as absorbing boundaries at the simulation volume borders to prevent unphysical back-scattering into the simulation volume [39]. Within a quazi-3D analysis, an ideal electromagnetic plane wave of circular polarization and angular momentum |m| = 1 is introduced into simulation volume using the TFSF approach [43]. To model focal spot shape due to scattering of a linearly-polarized plane wave (used in the experiment), we first performed two simulations using circularly-polarized plane waves of angular momenta m = 1 and m = −1. We then coherently added thus obtained solution fields to get a solution corresponding to a linearly polarized plane wave excitation field.
From Fig. 4, we notice an overall reasonable agreement between experimental and numerical data -namely, both reveal strongly sub-wavelength resolution of the THz SIL microscope beyond the Abbe diffraction limit in free space. The experimentally measured resolution parameter is found in the range of δ ≃ 0.15-0.4 for all measured samples that represent distinct modalities of the THz SIL operation as detailed earlier and summarized in Fig. 2(c). As follows from the earlier discussion, results of the 2D-FDTD simulations give the most optimistic lower bound of the SIL resolution parameter, which falls into the δ ≃ 0.15-0.2 range for all the samples. At the same time, quasi-3D-FEFD modeling predicts the values of a resolution parameter in the range of δ ≃ 0.2-0.4, which is bounded by the minimal and maximal sizes of an elliptical SIL focal spot.
We find that 3D simulations demonstrate inherently elliptical character of the THz SIL focused spot. In fact, beam spot is somewhat elongated in the direction of the electric field polarization E in the object plane. At the same time, we conclude that 2D and 3D analysis give almost identical results for the smaller size of a generally elliptical SIL focal spot, and, therefore, 2D simulations describe correctly the lower bound for the SIL sensitivity. As 2D simulations are generally somewhat faster than the quazi-3D ones, it might be a method of choice when estimating the best possible resolution of a SIL modality within a multi-parameter characterization/optimization studies. Results of the 2D simulations are also relevant for the experimental verification presented in this paper as planar metal/object interfaces for all samples stayed perpendicular to the semiminor axis of an elliptical focal spot during scanning, which should result in optimal resolution. At the same time, a larger size of an elliptical focal spot (not probed explicitly in our experiments) can only be predicted using 3D simulations, while also setting an upper bound on the SIL resolution. We, therefore, expect that in general, experimentally measured SIL resolution should fall between these two limits, and it should depend on the orientation and shape of the object boundary.
In fact, for the particular case of a flat object boundary positioned perpendicularly to the semi-small axis of the focal beam spot we expect resolution close to the lower bound of the SIL resolution range. At the same time, measured worse-than-expected resolution is probably due to misalignment of a sample boundary (should be strictly perpendicular to the focal spot major axis), as well as due to scattering at the boundary shape and material irregularities due to mechanical nature of the sample preparation. The question of resolution dependence on the object boundary shape is, however, beyond the scope of this paper and deserves a separate study. Although some differences are observed between the results of experimental studies and numerical simulations, the resolution of our THz SIL microscope appeared to be object dependent in all cases, and falling within the well understood limits set by the elliptical focal spot sizes. Moreover, for all considered materials, the SIL microscope provides strongly sub-wavelength resolution, thus, validating the use of this novel THz imaging modality in various applications.

Discussions
The detailed THz SIL microscope can be effectively used in a variety of applications, where sub-wavelength spatial resolution and high energy efficiency are required. Particularly, in pharmaceutical industry, THz SIL microscopy can be applied to study heterogeneity of chemical composition, porosity, and coatings of tablets [56,57]. In consumer electronics, it can be useful for non-destructive testing of electrical circuit boards [19,58]. Recently, THz spectroscopy and imaging was successfully applied to diagnosis of malignant and benign neoplasms with different nosologies and localizations [59]. Thanks to the high spatial resolution, the THz SIL microscopy is able to detect small-size pathologies and improve accuracy of the tumor margins delineation. In all these applications, the optical properties of imaging objects are expected to be within the considered range. Thus, the THz SIL microscopy is expected to provide high resolution in all these applications.
Due to linear nature of the Maxwell's equations [40], results presented in this work are general and can be transferred to analysis of SIL optical systems operating in other spectral ranges, such as visible and infrared. At the same time, thus discussed limitations may not be directly transferrable to other near-field imaging modalities that rely on different physical mechanisms; for example, see Refs. [28][29][30][31][32][33][34][35]. In order to uncover the object-dependent performance of other imaging modalities, one should perform independent analysis, involving analytical or computational methods of electrodynamics, while taking into account the particular physical modalities and hardware peculiarities responsible for the imaging system high spatial resolution.
Further studies of the object-dependent resolution of SIL microscopy remains an interesting research topic. While 3D simulations of the SIL performance involving different methods of computational electrodynamics seem to be a daunting task, they could indeed provide much better prediction for the spatial resolution of the SIL optical systems, as is evident from Fig. 4. At the same time, using efficient quazi-3D simulations is limited to rotationally symmetric systems, while further research in SIL resolution sensitivity to the object shape should be rather performed within a truly 3D simulation framework. We also note that the object inhomogeneous nature (fluctuations of the sample optical properties over its aperture) and related effects of the Mie scattering [41] of the evanescent and propagating waves can have an impact on the SIL microscopy resolution. Study of such phenomena, however, is out of the scope of the present paper and is also deferred to the future work.
The described SIL imaging modality should not be confused with the dielectric microparticleassisted microscopy [60][61][62], which usually relies on the photonic jet effect (i.e. the electromagnetic beam confinement behind a mesoscale dielectric particle, ≥ λ), or even the whispering gallery modes [40,[63][64][65][66]. To date, such imaging modalities have been extensively studied in the visible and infrared ranges, being capable of the competitive resolution of ∼ 0.25-0.5λ; while, recently, they were also transferred into the THz range [67,68]. The main disadvantage of these techniques is the necessity of using specialty particle patterned substrates at the imaging plane. Finally, as pointed out in Ref. [69], a favorable combination of the solid immersion and photonic jet phenomena is possible in high-refractive-index dielectric particles, which can further boost the THz microscopy performance. Such a combination forms a topic for separate full-blown study.

Conclusions
In this paper, numerical and experimental studies of the object-dependent spatial resolution of the original reflection-mode THz SIL microscopy system were carried out. It is revealed that the resolution of our system remains strongly sub-wavelength δ ≃ 0.15-0.4 for a wide range of sample refractive indices n obj = 1.0-5.0 and power absorption coefficients α obj = 0-400 cm −1 , thus, manifesting strong potential of the THz SIL microscopy for various applications.