## Abstract

Solid Immersion (SI) microscopy is a modern imaging modality that overcomes the Abbe diffraction limit and offers novel applications in various branches of visible, infrared, terahertz, and millimeter-wave optics. Despite the widespread use, SI microscopy usually results in qualitative imaging. Indeed, it presents only the raw distributions (in the image plane) of the backscattered field intensity, while unlocking the information about the physical properties of an imaged object, such as its complex refractive index (RI) distribution, requires resolving the inverse problem and remains a daunting task. In this paper, a method for resolving the SI microscopy inverse problem is developed, capable of reconstructing the RI distribution at the object imaging plane with subwavelength spatial resolution, while performing only intensity measurements. The sample RI is retrieved via minimization of the error function that characterizes discrepancy between the experimental data and the predictions of analytical model. This model incorporates all the key features of the electromagnetic-wave interaction with the SI lens and an imaged object, including contributions of the evanescent and ordinary-reflected waves, as well as effects of light polarization and wide beam aperture. The model is verified numerically, using the finite-element frequency-domain method, and experimentally, using the in-house reflection-mode continuous-wave terahertz SI microscope. Spatial distributions of the terahertz RIs of different low-absorbing optical materials and highly absorbing biological objects were studied and compared to *a priori* known data to demonstrate the potential of the novel SI microscopy modality. Given the linear nature of the Maxwell’s equations, the developed method can be applied for subwavelength-resolution SI microscopy at other spectral ranges.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Solid Immersion (SI) microscopy was developed in 1990 to overcome the Abbe diffraction limit in the visible range [1]. It exploits the SI effect, the essence of which is reduction in the dimensions of a focal spot formed in free space, a small distance away from the flat surface of a SI lens with a high refractive index (RI). Field contributions from both the transmitted (propagating) waves and evanescent waves, which originate due to the otal internal reflection (TIR) phenomenon, have to be taken into account when describing the SI lens focusing [2]. Generally, spatial resolution of the lens- or mirror-based imaging systems can be enhanced by a factor of the SI lens RI, while such lens features high optical throughout, thanks to the lack of any subwavelength apertures or probes in the optical path.

Such a superior performance has paved the way for the widespread use of SI microscopy in various branches of the visible and infrared optics. Particularly, it finds many applications in supper-resolution optical imaging [3,4], optical data storage [5], thermal emission microscopy [6], Raman imaging [7], photolitography [8], physics of semiconductors and superconductivity [9–12], non-destructive evaluations [13], biomedicine [14], single-photon emitters, and quantum optics [15]. Simultaneously, mechanically robust monolithic SI lenses and mirrors were developed with the help of various microfabrication techniques [16–21]. Recently, the principles of SI microscopy were translated into the millimeter-wave and terahertz (THz) spectral ranges [22–29]. Particularly, in the THz range (with the frequencies of $\nu \sim {10^{- 1}} {-} {10^1}\;{\rm THz}$), our team achieved the spatial resolution of $0.15\lambda$ ($\lambda$ is a free-space wavelength) for the reflection-mode continuous-wave (CW) SI microscopy system, as well as adapted it for imaging of amorphous objects and soft biological tissues [29–31].

Electromagnetic-wave focusing by different SI lens arrangements was studied theoretically in Refs. [32–35], where polarization and apodization of the electromagnetic beam, wavefront aberrations, interplay between propagating and evanescent waves, manufacturing errors, and misalignments in the optical path were carefully analyzed. Since SI lenses feature wide numerical apertures and operate in the near-field regime, one has to resort to the computational electrodynamics to accurately predict performance of such optical systems, including estimations of the focal spot geometry, spatial resolution, and depth of field. In Ref. [31], the effect of the sample complex RI on the SI lens performance was studied numerically, using a combination of the three-dimensional finite-element frequency-domain (3D FEFD) and the two-dimensional finite-difference time-domain methods, as well as experimentally, using a reflection-mode CW terahertz (THz) SI microscope. While the SI lens resolution was found to be object dependent and non-isotropic, it remained strongly subwavelength and below the Abbe diffraction limit for a majority of the practically important samples [31].

Despite its widespread use, SI microscopy usually results in qualitative imaging. Indeed, in most cases, SI images present only the raw distributions (in the image plane) of the backscattered field intensity. Meanwhile, unlocking information about the physical properties of an imaged object, such as its complex RI, remains a daunting task, since it requires resolving the SI microscopy inverse problem. Standard approaches for estimation of the sample optical properties, which are used in the dielectric spectroscopy, typically rely on the plane-wave, paraxial-beam, and far-field approximations [36]. They cannot be applied directly to the analysis of a SI lens due to its wide aperture and a subtle interplay between the evanescent and propagating waves in the near-field zone behind the SI lens flat surface. This pushes further research into resolving the inverse problem of SI microscopy and requires the fully vectorial solutions in the near field.

To address this challenge, in this paper, a novel method for resolving the inverse problem of SI microscopy is reported that aims at reconstructing the sample RI distribution with a subwavelength spatial resolution. In this method, the sample RI profile is found via minimization of an error functional that quantifies discrepancy between the experimental data and the analytical model. This model comprehensively describes the SI lens reflectivity as a function of the sample optical properties. It is first verified numerically, using the 3D FEFD modeling, and then experimentally, using the in-house reflection-mode CW THz SI microscope. The developed method is then used to study two practically important object types—i.e., low- and highly absorbing media. Finally, as a demonstration of utility of our methodology, RI distributions across the imaging plane are reconstructed for the freshly excised intact tissues and glioma model from rats *ex vivo*. The quantitative SI microscopy reported in this paper greatly expands the capabilities of this imaging modality in such demanding applications as condensed matter physics, material sciences, non-destructive testing, biology, and medicine.

## 2. RESULTS

#### A. Functioning of a THz SI Microscope

In this study, we consider an original THz SI lens from Refs. [29–31]. This lens is a key element of the in-house reflection-mode CW THz SI microscope [29–31], which is detailed in the supplemental document.

Our THz Si lens comprises a wide-aperture aspherical singlet (made of high-density polyethylene, HDPE) and a near-focal high-refractive-index hemispherical lens (made of high-resistivity float-zone silicon, HRFZ-Si). The latter serves as a resolution enhancer, allowing us to improve the resolution of basic HDPE aspherical singlet by a factor of HRFZ-Si RI at THz frequencies ${n_{{\rm Si}}} \simeq 3.415$ [29,37]. Moreover, HRFZ-Si possesses negligible absorption and chromatic dispersion in the THz range, which underlay the high energy efficiency of our SI lens as well as its potential for multispectral THz applications. The considered SI lens features the maximal aperture angle of ${\theta _{\max}} \simeq {40^ \circ}$.

Despite the HRFZ-Si hemisphere operating as a unitary optical element, it is comprised of the two mechanically separate components—namely, a rigidly fixed HRFZ-Si hypohemisphere and a movable HRFZ-Si window. As detailed in Ref. [29] and the supplemental document, such a composite construction of the HRFZ-Si hemisphere makes possible imaging of amorphous objects and soft biological tissues, thus significantly extending the capabilities of the THz SI microscopy in biophotonics.

This microscope is capable of imaging of both solid objects and soft biological tissues *ex vivo* with the spatial resolution as high as $r/\lambda \in ({0.15,0.4})$, which depends on the optical properties of an imaged media and appears to be asymmetric in the image plane but remains strongly subwavelength for a variety of practically important objects [31]. The described SI microscope also features small depth of field $\delta \simeq 0.1 {-} 0.2\lambda$ [30].

In our THz SI microscope, the sample is handled behind a movable HRFZ-Si window and imaged by a raster scan of its surface with a focused THz beam. For this aim, the HRFZ-Si window is mounted on a motorized 2D linear translation stage, which allows displacements in the lateral plane without affecting optical alignment [29,30]. In this way, a spatial distribution of the backscattered field intensity $I({\textbf r})$ is detected (${\textbf r}$ is a radius vector at the imaging plane), thus forming an initial condition for resolving the SI microscopy inverse problem.

#### B. General Approach to Resolve the SI Microscopy Inverse Problem

Consider an object with some spatial distribution of the complex RI across the image plane $\tilde n({\textbf r})$, which can be defined at each point ${\textbf r}$ as

where $n$ and $n^{\prime \prime} $ are its real and imaginary parts, $\alpha$ is the power absorption coefficient in [${{\rm cm}^{- 1}}$], and ${c_{0}} \simeq 3 \times {10^8}\;{\rm m}/{\rm s}$ is the speed of light in free space. Also, consider the following three assumptions regarding the thickness $l({\textbf r})$, heterogeneity, and roughness of an imaged object:- • At each point ${\textbf r}$, the object thickness $l$ is much larger than the depth of field $\delta$ of a SI lens $l \ge \delta \simeq 0.1\lambda$, which allows us to neglect scattering of waves from the back surface of an object, along with related standing waves inside an object.
- • Optical properties of an imaged object vary slowly over the image plane on a scale set by the SI lens resolution $r/\lambda \in ({0.15,0.4})$ and the depth of field $\delta \simeq 0.1\lambda$, which allows us to consider an object to be homogeneous within the focal spot.
- • Object surface roughness (in contact with the SI lens) is much smaller than the THz wavelength ($\ll \lambda$), the spatial resolution ($\ll r$), and the depth of field ($\ll \delta$), while an interface between the HRFZ-Si window and an object can be assumed to be perfectly flat. This way, roughness at the SI/object interface would not impact the performance of a SI system.

In other words, we consider the imaged object to be homogeneous within the THz beam spot and the depth of field of our THz SI microscope, which allows us to describe its THz response using the effective optical properties.

Under these assumptions, reconstruction of the object RI distribution $\tilde n({\textbf r})$ across the imaging plane can be achieved by minimizing an error function that quantifies a discrepancy between the experimental data and the theoretical model:

where ${H_{\exp}}$, ${H_{{\rm th}}}$ are the experimental and theoretical transfer functions. The function ${H_{\exp}}$ is defined based on the experimental data as where $I_{\exp}^{{\rm obj}}$ and $I_{\exp}^{{\rm ref}}$ are the intensities of the electromagnetic fields that are backscattered from the SI lens. $I_{\exp}^{{\rm obj}}$ is obtained for an imaged object, while $I_{\exp}^{{\rm ref}}$ is obtained for a reference medium with known optical properties. In this work, we use air (a SI lens without an object) as a reference medium with ${n_{{\rm ref}}} = 1.0$ and ${\alpha _{{\rm ref}}} = 0$. The function ${H_{{\rm th}}}$ is defined in a similar manner:An important issue for the successful resolution of the SI microscopy inverse problem is defining an intuitive theoretical model that describes the object-dependent backscattered signal ${I_{{\rm th}}}$, which, at the same time, accounts for all the principal physical mechanisms of the wave backscattering from the SI lens with an object at its shadow side. Such a model should be also valid for a wide range of practically important object RIs $n$ and losses $\alpha$. Since most biological tissues, liquids, polymers, composite structures, as well as majority of crystalline media have quite low RIs and finite absorption coefficients, in this work, we limited our analysis to the practically important ranges of $n \in ({1.0,5.0})$ and $\alpha \in ({{{0,10}^3}})\;{{\rm cm}^{- 1}}$.

Thanks to rather simple geometry of our SI lens, we resort to deriving the fully analytical and, thus, computationally efficient model, which accounts for such important effects as interplay between the ordinary reflected and TIR waves, light polarization, wide beam aperture, near-field operation regime, as well as standing waves inside the Si hemisphere.

#### C. Analytical Model of the SI Lens

In Fig. 1, we show schematic of the electromagnetic-wave reflection from the SI lens with some geometrical definitions used for deriving the analytical model of the reflected signal ${I_{{\rm th}}}({\tilde n})$. It is worth noting that the HDPE singlet is not considered explicitly, while it enters the formulation only via the numerical aperture defined by ${\theta _{\max}}$.

In our analytical model, the function ${I_{{\rm th}}}({\tilde n})$ is defined as

From Fig. 1, the ordinary reflected waves with $ s $ and $ p $ polarizations exist at incidence angle less than the object-dependent critical TIR angle $\theta \lt {\theta _{{\rm TIR}}}(n) = {\arcsin} ({n/{n_{{\rm Si}}}})$, while oppositely, the TIR reflectance occurs at the incidence angles $\theta \ge {\theta _{{\rm TIR}}}$. For $ s $ and $ p $ polarizations, the corresponding net amplitude of the reflected wave is given by

For the infinite coherence length of our CW THz emitter, the reflectivity ${R^{{{s/p}}}}$ is defined as

General formulations [Eqs. (2)–(4)], aided by the described analytical model [Eqs. (5)–(12)], define a relation between the measured images and the object RI distribution and can be applied to resolve the SI microscopy inverse problem.

While the analytical model ${I_{{\rm th}}}$ offers clear physical understanding of the Si lens reflectivity properties, its accuracy might be somehow limited due to a number of approximations used in its derivation. Particularly, our model is formulated for an infinite object thickness. It also assumes both ordinary reflection and TIR at the intersections of an optical axis with the lens interface. Finally, the aspherical singlet enters our analytical analysis only via its numerical aperture ${\theta _{\max}}$, while the beam attenuation and apodization, that might occur due to the Fresnel losses, variable thickness, and THz-wave absorption in the singlet are neglected.

In the analytical and numerical modeling of SI lens reflectivity, both amplitude and phase of the electromagnetic field as well as lossy objects are considered. However, since most SI microscopes only detect the field intensity, an important practical question is whether it is possible to resolve the inverse problem for the object complex RI [Eqs. (2)–(4)] using only amplitude data. Clearly, using only amplitude information somewhat limits the generality of the developed method, since it is impossible to simultaneously estimate the two independent parameters $n$, $\alpha$ based on a single measured value ${I_{\exp}}$. At the same time, in what follows we show how amplitude data alone can be sufficient to recover the real part of the object RI in several cases, such as low-loss samples, or lossy objects, where the relation between the real and imaginary parts of their RI can be reliably modeled.

#### D. Numerical Modeling of the SI Lens

In order to verify the analytical model, we use the 3D FEFD method as detailed in the supplemental document. It allows us to numerically model the object-dependent backscattered field intensity ${I_{{\rm th}}}$ and then obtain numerical predictions for the theoretical transfer function ${H_{{\rm th}}}$, while considering both infinitely thick ($l \to \infty$) lossless ($\alpha = 0$) objects and those with finite thickness $l$ and loss $\alpha$. Computational-intensive full-vectorial numerical analysis relies on the first principles and thus accounts for all peculiarities of the electromagnetic-wave–SI lens interaction.

In Fig. 2(a), analytical and numerical functions ${H_{{\rm th}}}$ are compared for the technical parameters of our SI lens, a variable $n$, zero loss $\alpha = 0$, and infinite thickness $l \to \infty$ of an object. The developed method should provide a unique inverse problem solution—i.e., ${H_{{\rm th}}}$ should determine unambiguously the relation between an image intensity ${I_{\exp}}$ and an object RI $n$. This demands from ${H_{{\rm th}}}$ to have a monotonic behavior in the desired range of RI variations. In Fig. 2(b), we show the first derivative $d{H_{{\rm th}}}/dn$ of the analytical and numerical functions. The analytical curve possesses a negative first derivative and thus decreases monotonically in the RI range of $n \in ({1.0,4.5})$, thus providing a unique inverse problem solution. The numerical data overall agree with the analytical estimates and also behave largely monotonically. At lower object RIs $n \le 2.2$, the TIR effect has significant contribution at the Si/object interface, the evanescent waves are excited [31], and pronounced changes in ${H_{{\rm th}}}$ are observed, thus leading to higher sensitivity of the developed method to small RI variations. At higher RIs $n \gt 2.2$, the TIR condition at the interface is violated, and only the ordinary propagating waves exist [31], thus leading to the less-pronounced changes in ${H_{{\rm th}}}$ and smaller sensitivity of the RI reconstruction.

In Figs. 2(c) and 2(d), the analytical transfer function ${H_{{\rm th}}}$ is compared with the numerical one for the infinitely thick ($l \to \infty$) object with RIs of $n \in ({1.0,5.0})$ and power absorption coefficient of $\alpha \in ({{{10}^{- 1}}{{,10}^3}})\;{{\rm cm}^{- 1}}$. Here, the analytical model shows good quantitative agreement with the numerical data. The observed discrepancy between analytical and numerical data can be attributed to a number of factors mostly related to various approximations used in the formulation of the analytical model. Also, in the 3D FEFD modeling, the backscattered field intensity is determined at a small distance in front of the HDPE singlet (see the supplemental document), while in the analytical model and experiments, the intensity is determined at the far-field limit. The overall tendency of reduced reflection when increasing object RI $n$ is related to the decrease of a critical angle of incidence ${\theta _{{\rm TIR}}}$ at the Si/object interface (see the supplementary document), resulting in higher transmission intensity through the Si hemisphere. Furthermore, the dielectric contrast at the interface decreases while $n$ increases up to ${n_{{\rm Si}}}$, which also leads to a SI lens reflectivity drop. Oppositely, a minor increase in ${H_{{\rm th}}}$ at high RIs $n \gt 4.5$ might be attributed to an enhanced dielectric contrast and thus higher Fresnel reflectivity of a system [31].

Finite thickness $l$ of an image object might complicate reconstruction of the RI profile as detailed in the supplemental document, where the transfer function ${H_{{\rm th}}}$ is calculated numerically using the 3D FEFD method for the object thicknesses of $l = 0.1\lambda$, $\lambda$, $10\lambda$, and $\infty$. When dealing with a lossless object ($\alpha = 0$), finite $l$ promotes standing waves inside an object and thus results in modulations in ${H_{{\rm th}}}$ that are more pronounced at higher RIs. In fact, due to the nature of standing waves, periodic variation in ${H_{{\rm th}}}$ as a function of $n$ is expected to have a characteristic period of $\Delta n \sim \lambda /2l$. While $l$ reduces, the modulations occupy wider RI range and limit the ranges of $n$ and $\alpha$, where the inverse problem can be unambiguously resolved. Finite $\alpha$ can attenuate these standing waves and suppress modulations in ${H_{{\rm th}}}$, when

As ${H_{{\rm th}}}$ is monotonic and does not change considerably as a function of loss when $\alpha \le 10\;{{\rm cm}^{- 1}}$ [Fig. 2(c)], the developed method based on the intensity-only measurements can be used for unambiguous estimation of the object RI $n$ for both lossless and low-loss objects, while no *a priori* knowledge about $\alpha$ is required. Moreover, for higher losses $\alpha \gt 10\;{{\rm cm}^{- 1}}$, due to smooth and largely monotonic behavior of ${H_{{\rm th}}}$, the same method can be applied to study complex RI $\tilde n$ of lossy media, assuming that some physically reasonable relation between $n$ and $\alpha$ can be defined. For example, in the case of lossy objects of known origin (such as mixes of liquids with different concentrations), the object RI and its loss are confined to a certain curve in the ($n$, $\alpha$) space, defining a specific pattern of the response function ${H_{{\rm th}}}$, that can be used for resolving the inverse problem.

#### E. Experimental Results

The developed method of the object RI estimation was then verified experimentally using intensity data of the in-house reflection-mode CW THz SI microscope, operating at $\nu = 0.6\; {\rm THz}$ or $\lambda \simeq 500\; \unicode{x00B5}{\rm m}$; see the supplemental document. In what follows, we separately consider imaging of low- and high-absorbing objects.

### 1. Low-Absorbing Media

First, several thick ($\gg \lambda$) slabs made of the materials with distinct *a priori* known RIs $n$ and relatively low absorption coefficients $\alpha \le 10\;{{\rm cm}^{- 1}}$ are analyzed. Among those are polymethylpenthene (TPX), HDPE, polymethyl methacrylate (PMMA), and crystalline quartz (${{\rm SiO}_{2}}$). During measurements, polished slabs were placed in direct contact with the Si window of a SI lens, without an air gap between them.

In Figs. 3(a)–3(h), results of the THz SI microscopy of these slabs are shown in the form of the THz intensity images $I({\textbf r})$ and the estimated RI distributions $n({\textbf r})$, while each THz image captures a fragment of a slab bordering air. RI distributions $n({\textbf r})$ over the imaged area are used to calculate the mean RI value and its standard deviation for each slab. In Fig. 3(i), we compare the found averaged parameters with their known values and find an excellent agreement.

We thus conclude that inversion algorithm for RI reconstruction based on transfer function ${H_{{\rm th}}}$ works well for low-loss materials and requires no *a priori* knowledge of the object loss $\alpha$. Thereby, the developed method can be directly applied for studying relatively low-loss objects in such demanding applications as THz material sciences, non-destructive testing of dielectric media, as well as chemistry and pharmacology.

#### F. Highly Absorbing Media

Next, we consider highly absorbing media, for which some *a priori* knowledge of absorption coefficient $\alpha$ or a physical relation between the RI $n$ and loss parameter $\alpha$ should be used.

Measurements of highly absorbing media are of particular importance in THz biophotonics, where the high dipole moment of ${\rm H}_2{\rm O}$ molecules plays a dominant role in the relaxation dynamics of the THz dielectric response of biological tissues and solutions of biomolecules [40]. High water content leads to strong THz-wave absorption by hydrated biological objects and thus limits the depth of THz-wave penetration in them. At the same time, tissue water is a main endogenous marker of pathological processes in tissues in label-free THz medical diagnosis [40]. For these reasons, in our study of highly absorbing media, we concentrate on measuring water content in biological tissues and aqueous mixtures. Particularly, by determining the spectral distribution of the tissue refractive index $n({\textbf r})$ using THz SI microscopy, we can then deduce the spatial distribution of local tissue water content $C({\textbf r})$ for further tissue characterization.

In our analysis, we use Ref. [41], where a model was introduced to describe the THz complex RI $\tilde n$ of tissues as a function of the tissue water content $C$ (per unitary volume). This model treats a tissue as a two-component media, comprised of tissue water, with its complex RI ${\tilde n_{{{\rm H}2{\rm O}}}}$, and all other less-polar tissue components (such as cellular matrix, cell components, proteins, lipids, deoxyribonucleic acid, amino acids, and a variety of other biological molecules) that feature the cumulative complex RI ${\tilde n_{{\rm dry}}}$, which corresponds to the dehydrated (dried or lyophilized) tissues. In place of a simplified linear decomposition approach, which was used to determine the tissue response in Ref. [41], in this study, we resort to the more physically rigorous Bruggeman model of the two-component media [42]:

In Fig. 4(a), we present the analytical transfer function ${H_{{\rm th}}}$ in the RI range of $n \in ({1.6,2.4})$ and the absorption coefficient range of $\alpha \in ({0,200})$. Using the measured complex RI of liquid water ${\tilde n_{{{\rm H}2{\rm O}}}} = 2.32 - i0.69$ (${\alpha _{{{\rm H}2{\rm O}}}} = 174\;{{\rm cm}^{- 1}}$) from the supplemental document and that of dehydrated brain tissues ${\tilde n_{{\rm dry}}} = 1.8 - i0.05$ (${\alpha _{{\rm dry}}} = 12\;{{\rm cm}^{- 1}}$) from Refs. [41,46], at 0.6 THz, in Fig. 4(a), as a white-color curve, we show an evolution of the brain tissue optical properties, with water content increasing from $C = 0$ to 1 as predicted by Eq. (14). One can notice the close-to-linear character of the theoretical dependence. Furthermore, in Fig. 4(a), we show optical properties at 0.6 THz of some representative healthy tissues and neoplasms from humans, which appear to be very close to the defined model of the tissue complex RI.

In what follows, we use Eq. (14) within the inverse problem formulation in order to relate the object RI $n$ and losses $\alpha$ within the context of characterization of the water content in biological objects. More precisely, using Eq. (14) for the water concentrations of $C \in ({0,1})$, from Fig. 4(a), we can then calculate a cross section of the analytical transfer function $H_{{\rm th}}^{{\rm tissue}}$ as shown in Fig. 4(b). We notice that $H_{{\rm th}}^{{\rm tissue}}$ is monotonic within the considered range of tissue RIs. Thus, the developed model can be applied for the unambiguous assessment of spatial distributions of the tissues optical properties $n({\textbf r})$, $\alpha ({\textbf r})$, and water content $C({\textbf r})$.

As a test system to verify the developed method, we used propylene glycol (PG) and its aqueous solutions with the concentrations of 10%–90% and the step of 10% (by volume), which forms a very approximate tissue-mimicking phantom. Optical constants of PG and its aqueous solutions were measured using the transmission-mode THz pulse spectrometer; see the supplemental document. From Fig. 4(a), we notice that pure PG possesses quite small RI $n$ and loss $\alpha$, which are similar to those of dehydrated tissues. Optical constants $n$ and $\alpha$ of the PG aqueous solutions increase almost linearly with the water content $C$, thus mimicking (to some extent) the tissue optical properties. For studying liquids with the reflection-mode CW THz SI microscope, a sample cuvette is placed atop of the SI lens [see the inset of Fig. 4(b)]. The inferred RIs $n$ of PG and its aqueous solutions are then compared with the values obtained by spectroscopic measurements [see Fig. 4(c)]. We observe excellent agreement between the data sets measured by the THz SI microscopy and THz pulsed spectroscopy, thus justifying the validity of the developed approach. By studying a wide variety of PG aqueous solutions of different concentrations, we demonstrated that the developed method should work well with both highly hydrated biological samples (such as brain tissues with the water content $\simeq 70\%$ [45]) and less-hydrated tissues (such as bones with the water content of $\simeq 30\%$ [47]).

Next, to demonstrate practical utility of the developed method, it is applied for imaging of the freshly excised intact brain tissues and glioma model 101.8 from rats *ex vivo* [48] as detailed in the supplemental document. In Fig. 5, the results of the quantitative THz SI microscopy of a healthy rat brain and a tumor model are compared. The estimated distributions of tissue RI $n({\textbf r})$ and absorption coefficient $\alpha ({\textbf r})$ agree with the average values from previous studies [41,45], which involved THz pulsed spectroscopy and did not differentiate between the white matter and gray matter, as well as other neurovascular structures due to the diffraction-limited resolution. In turn, earlier estimations of water content in brain tissues differ considerably, depending on the applied physical principle and the experimental technique, and usually appear in the range of 55%–85% [45]. Thus, as compared to other experimental techniques, our method gives overall reasonable estimates of $C({\textbf r})$ but might somewhat underestimate this parameter due to a simplified model of the tissue complex RI used in the interpretation of the SIL data. Thanks to the superior resolution of our THz SI microscope, considerable heterogeneity of the brain tissues at the THz-wavelength scale is clear from the THz images.

From Figs. 5(b)–5(e), a pronounced difference is observed between the white matter and gray matter (cortex), which is consistent with observations in Refs. [49,50]. Particularly, white matter possesses higher RIs $n$ and absorption coefficients $\alpha$, which can be attributed to increased content of tissue water $C$ (in the form of an electrolyte solution inside axons) and myelin (a lipid-rich substance surrounding the axons). In Figs. 5(g)–5(i), we also notice that a tumor possess higher RIs $n$ and absorption coefficient $\alpha$ than those of intact tissues, which is mostly due to the higher water content $C$ in a tumor [41,45]. Also, a considerable mesoscale heterogeneity of a tumor is notable, due to various factors, such as the tumor cells’ accumulations, vessels, necrotic debris, and hemorrhages.

#### G. Discussions

On the one hand, the observed results of the quantitative THz SI microscopy of the rat brain justify a label-free contrast between intact tissues and a tumor, which was earlier pointed by several research groups [40,41,45,48]. On the other hand, our THz images highlight considerable heterogeneity of brain tissues at the scale posed by the THz wavelengths. This poses important problems for studying the THz-wave scattering effects, the interplay between Rayleigh and Mie scattering in the THz range, and the impact of absorption and scattering processes on the THz-wave extinction in tissues, as well as for adapting the radiative transfer theory for the needs of THz biophotonics, which are still to be addressed [40,48,51,52].

Indeed, tissue heterogeneities and related THz-wave scattering effects might either complicate THz diagnosis of tumors (being a reason of strong variability of the tissues response in the THz range) or become a source of additional useful information for the tissue characterization and differentiation. For example, in Ref. [45], THz dielectric spectroscopy of intact tissues and human brain gliomas of the WHO Grades I–IV [notice that low-grade gliomas (I, II) are benign, while high-grade ones (III, IV) are malignant] did not reveal any useful information for the discrimination between gliomas of the different grades, including benign and malignant. In turn, such structural heterogeneities of a tumor, as areas of necrosis and hemorrhages, are inherent specifically to the WHO Grade IV glioma (i.e., glioblastoma), which is mimicked by the consider glioma model 101.8 [48,53]. Since such heterogeneities were visualized by the SI microscopy in this study, it is possible to use them as a source of additional information for the differentiation between glioblastoma and other benign and malignant gliomas. We hope than in our future work we could find other structural features of benign and malignant neoplasms that can be clearly visualized by our THz SI microscope and that can open the road toward their THz diagnosis.

It is worth noting that our THz SI system, aided by the developed method of RI estimation, allows us to visualize practically important mesoscale ($\sim{10^1} {-} {10^{- 1}}\lambda$) heterogeneities in the biological tissues. Such heterogeneities can also lead to Mie scattering of the THz waves in tissues, which is accompanied by strong variability of the scattering parameters of tissues depending on their composition, local optical properties, and morphology [54]. Thus, the developed tool can be quite effective when solving the problem of THz-wave scattering in tissues, as well as when adapting the Mie scattering theory and the radiative transfer theory for the THz biophotonics. At the same time, tissue heterogeneities with smaller dimensions ($\le {10^{- 1}}\lambda$) that cannot be resolved by our THz SI system should manifest purely Rayleigh scattering regime, with its isotropic scattering indicatrix. For such tissues, the effective THz dielectric response (averaged within the THz beam spot) can still be used to describe the THz wave–tissue interactions within the frameworks of the effective medium theory [40].

The developed quantitative super-resolution THz SI microscopy holds strong potential in a number of demanding fields in material sciences, non-destructive measurements, quality control, security, biology, chemical, and pharmaceutical industries. In all these applications, high spatial resolution and high energy efficiency (thanks to the absence of any subwavelength probes in the optical scheme), along with its ability to extract quantitative information about an imaged object (in the form of the THz RI $n$ and material loss $\alpha$, water content $C$, or other related physical quantities) represent a favorable combination of technical advantages of the developed technique over other modalities of THz spectroscopy, imaging, and sensing.

In our opinion, among all the potential applications of this method, THz biophotonics is particularly attractive due to the strong sensitivity of THz waves to the content and state of tissue water [40]. In fact, tissue water is a main endogenous marker of pathological processes in the label-free diagnosis of malignant and benign neoplasms with different nosologies and localization, viability and traumatic injuries of tissues, diabetes mellitus, and so on. The developed THz SI microscopy modality puts a strong emphasis on the quantitative diagnosis of tissues, where local tissue RI $n$ and water content $C$, can play the role of physically reasonable signs for tissue diagnosis [40].

We note that the developed method can be further generalized for direct imaging of the complex (absorbing) media using the field amplitude detectors. Particularly, by detecting both the intensity and phase of the reflected light in a SI microscope and then interpreting their values using the developed SI lens reflectivity model, one could unambiguously reconstruct the material complex RI distribution without any *a priori* assumptions on the imaged material properties. In practice, this can be realized, for example, by adapting the principles of THz pulsed imaging, photomixing, and heterodyne or homodyne detection for simultaneous measurements of the THz signal amplitude and phase.

Thanks to linearity of the Maxwell’s equations, the developed method can be translated to other spectral ranges, such as visible, near- and mid-infrared, where many SI microscopy systems have been demonstrated. It can also be adapted to other SI lens designs, as long as the reflectivity of such lenses could be modeled analytically or numerically. Naturally, one has to be conscious of spectrally dependent material properties of various imaging elements, such as lenses, mirrors, and beam splitters used in the SI microscope. Similarly, frequency specific properties of the electromagnetic-wave emitters, such as monochromaticity, coherence length, polarization state, beam quality, and so on, as well as the type of the detected signal, such as field amplitude or power flux, should be accurately taken into account when developing the SI microscope or constructing the SI lens reflectivity model. Finally, it is worth noting that different electro- and magneto-dipole excitations of matter underlie the linear and nonlinear electrodynamic response of an imaged object in different spectral ranges. Among them, only linear dielectric polarization of matter and related phenomena can be studied using the developed method. Such a variability of physical processes can lead to considerable differences in the optical properties of an object, when it is imaged in different ranges of the electromagnetic spectrum. Such variations of the sample optical properties should be carefully taken into account when adapting the developed method to SI systems from other spectral ranges.

## 3. CONCLUSIONS

In this work, a method for reconstruction with subwavelength resolution of the spatial distribution of the material complex RI within the SI microscopy paradigm was developed and studied analytically, numerically, and experimentally. By measuring several test optical material combinations, as well as complex biological tissues *ex vivo*, using the in-house reflection-mode CW THz SI microscope, the predictive potential of the developed method was objectively verified. Given the linear nature of the Maxwell’s equations, the developed approach can be directly applied for subwavelength-resolution SI microscopy at other spectral ranges.

## Funding

Russian Science Foundation (RSF) (17-79-20346); Ministry of Science and Higher Education of the Russian Federation (075-15-2020-790); Canada Research Chairs (Ubiquitous Terahertz Photonics Project).

## Acknowledgment

Development of a novel method for solving the inverse problem of the SI microscopy and its experimental implementations for imaging of solid state objects, by K.I. Zaytsev and N.V. Chernomyrdin were supported by the Russian Science Foundation (RSF), Project # 17-79-20346. Fabrication of test objects, and the work with biological tissues and liquids, by G.M. Katyba, were supported by the Ministry of Science and Higher Education of the Russian Federation, Project # 075-15-2020-790. Numerical analysis of the THz SI lens performance using the FEFD method, by M. Skorobogatiy, was supported by the Canada Research Chairs Program, the Ubiquitous Terahertz Photonics Project. The authors are thankful to Pavel V. Nikitin from Sechenov University (Russia) and Anna I. Alekseeva from the Research Institute of Human Morphology (Russia) for their help with the tissue preparation and histology. They are also thankful to Dr. Irina N. Dolganova and Prof. Vladimir N. Kurlov from the ISSP RAS for the fruitful discussions.

## Disclosures

The authors declare no conflict of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## Supplemental document

See Supplement 1 for supporting content.

## REFERENCES

**1. **S. Mansfield and G. Kino, “Solid immersion microscope,” Appl. Phys. Lett. **57**, 2615–2616 (1990). [CrossRef]

**2. **T. Milster, J. Jo, and K. Hirota, “Roles of propagating and evanescent waves in solid immersion lens systems,” Appl. Opt. **38**, 5046–5057 (1999). [CrossRef]

**3. **D. Kang, C. Pang, S. M. Kim, H. S. Cho, H. S. Um, Y. W. Choi, and K. Suh, “Shape-controllable microlens arrays via direct transfer of photocurable polymer droplets,” Adv. Mater. **24**, 1709–1715 (2012). [CrossRef]

**4. **W. Fan, B. Yan, Z. Wang, and L. Wu, “Three-dimensional all-dielectric metamaterial solid immersion lens for subwavelength imaging at visible frequencies,” Sci. Adv. **2**, e1600901 (2016). [CrossRef]

**5. **B. Terris, H. Mamin, D. Rugar, W. Studenmund, and G. Kino, “Near-field optical data storage using a solid immersion lens,” Appl. Phys. Lett. **65**, 388–390 (1994). [CrossRef]

**6. **S. Ippolito, S. Thorne, M. Eraslan, B. Goldberg, M. Unlu, and Y. Leblebici, “High spatial resolution subsurface thermal emission microscopy,” Appl. Phys. Lett. **84**, 4529–4531 (2004). [CrossRef]

**7. **G. Lerman, A. Israel, and A. Lewis, “Applying solid immersion near-field optics to Raman analysis of strained silicon thin films,” Appl. Phys. Lett. **89**, 223122 (2006). [CrossRef]

**8. **L. Ghislain, V. Elings, K. Crozier, S. Manalis, S. Minne, K. Wilder, G. Kino, and C. Quate, “Near-field photolithography with a solid immersion lens,” Appl. Phys. Lett. **74**, 501–503 (1999). [CrossRef]

**9. **Q. Wu, R. Grober, D. Gammon, and D. Katzer, “Imaging spectroscopy of two-dimensional excitons in a narrow GaAs/AlGaAs quantum well,” Phys. Rev. Lett. **83**, 2652–2655 (1999). [CrossRef]

**10. **Z. Liu, B. Goldberg, S. Ippolito, A. Vamivakas, M. Unlu, and R. Mirin, “High resolution, high collection efficiency in numerical aperture increasing lens microscopy of individual quantum dots,” Appl. Phys. Lett. **87**, 071905 (2005). [CrossRef]

**11. **M. Sartison, S. Portalupi, T. Gissibl, M. Jetter, H. Giessen, and P. Michler, “Combining in-situ lithography with 3D printed solid immersion lenses for single quantum dot spectroscopy,” Sci. Rep. **7**, 39916 (2017). [CrossRef]

**12. **A. Bogucki, L. Zinkiewicz, M. Grzeszczyk, W. Pacuski, K. Nogajewski, T. Kazimierczuk, A. Rodek, J. Suffczyński, K. Watanabe, T. Taniguchi, P. Wasylczyk, M. Potemski, and P. Kossacki, “Ultra-long-working-distance spectroscopy of single nanostructures with aspherical solid immersion microlenses,” Light Sci. Appl. **9**, 48 (2020). [CrossRef]

**13. **S. Ippolito, P. Song, D. Miles, and J. Sylvestri, “Angular spectrum tailoring in solid immersion microscopy for circuit analysis,” Appl. Phys. Lett. **92**, 101109 (2008). [CrossRef]

**14. **L. Wang, B. Bateman, L. Zanetti-Domingues, A. Moores, S. Astbury, C. Spindloe, M. Darrow, M. Romano, S. Needham, K. Beis, D. Rolfe, D. Clarke, and M. Martin-Fernandez, “Solid immersion microscopy images cells under cryogenic conditions with 12 nm resolution,” Commun. Biol. **2**, 74 (2019). [CrossRef]

**15. **T. Schroder, F. Gadeke, M. Banholzer, and O. Benson, “Ultrabright and efficient single-photon generation based on nitrogen-vacancy centres in nanodiamonds on a solid immersion lens,” New J. Phys. **13**, 055017 (2011). [CrossRef]

**16. **D. Fletcher, K. Crozier, C. Quate, G. Kino, K. Goodson, D. Simanovskii, and D. Palanker, “Near-field infrared imaging with a microfabricated solid immersion lens,” Appl. Phys. Lett. **77**, 2109–2111 (2000). [CrossRef]

**17. **R. Brunner, M. Burkhardt, A. Pesch, O. Sandfuchs, M. Ferstl, S. Hohng, and J. White, “Diffraction-based solid immersion lens,” J. Opt. Soc. Am. A **21**, 1186–1191 (2004). [CrossRef]

**18. **I. Golub, “Solid immersion axicon: maximizing nondiffracting or Bessel beam resolution,” Opt. Lett. **32**, 2161–2163 (2007). [CrossRef]

**19. **D. Mason, M. Jouravlev, and K. Kim, “Enhanced resolution beyond the Abbe diffraction limit with wavelength-scale solid immersion lenses,” Opt. Lett. **35**, 2007–2009 (2010). [CrossRef]

**20. **M.-S. Kim, T. Scharf, M. Haq, W. Nakagawa, and H. Herzig, “Subwavelength-size solid immersion lens,” Opt. Lett. **36**, 3930–3932 (2011). [CrossRef]

**21. **T.-Y. Huang, R. Grote, S. Mann, D. Hopper, A. Exarhos, G. Lopez, G. Kaighn, E. Garnett, and L. Bassett, “A monolithic immersion metalens for imaging solid-state quantum emitters,” Nat. Commun. **10**, 2392 (2019). [CrossRef]

**22. **A. Pimenov and A. Loidl, “Focusing of millimeter-wave radiation beyond the Abbe barrier,” Appl. Phys. Lett. **83**, 4122–4124 (2003). [CrossRef]

**23. **B. Gompf, M. Gerull, T. Müller, and M. Dressel, “THz-micro-spectroscopy with backward-wave oscillators,” Infrared Phys. Technol. **49**, 128–132 (2006). [CrossRef]

**24. **B. Gompf, M. Dressel, H. Heer, and S. Martens, “THz-near-field micro-spectroscopy with backward-wave oscillators and a photo-induced aperture,” in *Joint 31st International Conference on Infrared Millimeter Waves and 14th International Conference on Terahertz Electronics* (2006), p. 17.

**25. **S. Kim, H. Murakami, and M. Tonouchi, “Transmission-type laser THz emission microscope using a solid immersion lens,” IEEE J. Sel. Top. Quantum Electron. **14**, 498–504 (2008). [CrossRef]

**26. **N. Chernomyrdin, A. Schadko, S. Lebedev, V. Tolstoguzov, V. Kurlov, I. Reshetov, I. Spektor, M. Skorobogatiy, S. Yurchenko, and K. Zaytsev, “Solid immersion terahertz imaging with sub-wavelength resolution,” Appl. Phys. Lett. **110**, 221109 (2017). [CrossRef]

**27. **J. Choi, *High-Speed Devices and Circuits with THz Applications*, 1st ed. (CRC Press, 2017).

**28. **A. Jomaso, C. Palacios, J. Mejía, R. Perez, and N. Qureshi, “Continuous wave microscopy based on solid immersion lens,” in *44th International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz)* (2019), pp. 1–2.

**29. **N. Chernomyrdin, A. Kucheryavenko, G. Kolontaeva, G. Katyba, I. Dolganova, P. Karalkin, D. Ponomarev, V. Kurlov, I. Reshetov, M. Skorobogatiy, V. Tuchin, and K. Zaytsev, “Reflection-mode continuous-wave 0.15λ-resolution terahertz solid immersion microscopy of soft biological tissues,” Appl. Phys. Lett. **113**, 111102 (2018). [CrossRef]

**30. **N. Chernomyrdin, V. Zhelnov, A. Kucheryavenko, I. Dolganova, G. Katyba, V. Karasik, I. Reshetov, and K. Zaytsev, “Numerical analysis and experimental study of terahertz solid immersion microscopy,” Opt. Eng. **59**, 061605 (2019). [CrossRef]

**31. **V. Zhelnov, K. Zaytsev, A. Kucheryavenko, G. Katyba, I. Dolganova, D. Ponomarev, V. Kurlov, M. Skorobogatiy, and N. Chernomyrdin, “Object-dependent spatial resolution of the reflection-mode terahertz solid immersion microscopy,” Opt. Express **29**, 3553–3566 (2021). [CrossRef]

**32. **M. Baba, T. Sasaki, M. Yoshita, and H. Akiyama, “Aberrations and allowances for errors in a hemisphere solid immersion lens for submicron-resolution photoluminescence microscopy,” J. Appl. Phys. **85**, 6923–6925 (1999). [CrossRef]

**33. **K. Karrai, X. Lorenz, and L. Novotny, “Enhanced reflectivity contrast in confocal solid immersion lens microscopy,” Appl. Phys. Lett. **77**, 3459–3461 (2000). [CrossRef]

**34. **L. Helseth, “Roles of polarization, phase and amplitude in solid immersion lens systems,” Opt. Commun. **191**, 161–172 (2001). [CrossRef]

**35. **C. Liu and S.-H. Park, “Numerical analysis of an annular-aperture solid immersion lens,” Opt. Lett. **29**, 1742–1744 (2004). [CrossRef]

**36. **I. Pupeza, R. Wilk, and M. Koch, “Highly accurate optical material parameter determination with THz time-domain spectroscopy,” Opt. Express **15**, 4335–4350 (2007). [CrossRef]

**37. **D. Grischkowsky, S. Keiding, M. van Exter, and C. Fattinger, “Far-infrared time-domain spectroscopy with terahertz beams of dielectrics and semiconductors,” J. Opt. Soc. Am. B **7**, 2006–2015 (1990). [CrossRef]

**38. **W. Withayachumnankul, B. Fischer, and D. Abbott, “Material thickness optimization for transmission-mode terahertz time-domain spectroscopy,” Opt. Express **16**, 7382–7396 (2008). [CrossRef]

**39. **M. Zhai, A. Locquet, and D. Citrin, “Pulsed THz imaging for thickness characterization of plastic sheets,” NDT&E Int. **116**, 102338 (2020). [CrossRef]

**40. **K. Zaytsev, I. Dolganova, N. Chernomyrdin, G. Katyba, A. Gavdush, O. Cherkasova, G. Komandin, M. Shchedrina, A. Khodan, D. Ponomarev, I. Reshetov, V. Karasik, M. Skorobogatiy, V. Kurlov, and V. Tuchin, “The progress and perspectives of terahertz technology for diagnosis of neoplasms: a review,” J. Opt. **22**, 013001 (2019). [CrossRef]

**41. **S. Yamaguchi, Y. Fukushi, O. Kubota, T. Itsuji, T. Ouchi, and S. Yamamoto, “Origin and quantification of differences between normal and tumor tissues observed by terahertz spectroscopy,” Phys. Med. Biol. **61**, 6808–6820 (2016). [CrossRef]

**42. **V. Markel, “Introduction to the Maxwell Garnett approximation: tutorial,” J. Opt. Soc. Am. A **33**, 1244–1256 (2016). [CrossRef]

**43. **E. Pickwell, A. Fitzgerald, B. Cole, P. Taday, R. Pye, T. Ha, M. Pepper, and V. Wallace, “Simulating the response of terahertz radiation to basal cell carcinoma using ex vivo spectroscopy measurements,” J. Biomed. Opt. **10**, 064021 (2005). [CrossRef]

**44. **P. Ashworth, E. Pickwell-MacPherson, E. Provenzano, S. Pinder, A. Purushotham, M. Pepper, and V. Wallace, “Terahertz pulsed spectroscopy of freshly excised human breast cancer,” Opt. Express **17**, 12444–12454 (2009). [CrossRef]

**45. **A. Gavdush, N. Chernomyrdin, G. Komandin, I. Dolganova, P. Nikitin, G. Musina, G. Katyba, A. Kucheryavenko, I. Reshetov, A. Potapov, V. Tuchin, and K. Zaytsev, “Terahertz dielectric spectroscopy of human brain gliomas and intact tissues ex vivo: double-Debye and double-overdamped-oscillator models of dielectric response,” Biomed. Opt. Express **12**, 69–83 (2021). [CrossRef]

**46. **K. Lee, K. Jeoung, S. Kim, Y.-B. Ji, H. Son, Y. Choi, Y.-M. Huh, J.-S. Suh, and S. Oh, “Measuring water contents in animal organ tissues using terahertz spectroscopic imaging,” Biomed. Opt. Express **9**, 1582–1589 (2018). [CrossRef]

**47. **S. Freer, C. Sui, S. Hanham, L. Grover, and M. Navarro-Cía, “Hybrid reflection retrieval method for terahertz dielectric imaging of human bone,” Biomed. Opt. Express **12**, 4807–4820 (2021). [CrossRef]

**48. **A. Kucheryavenko, N. Chernomyrdin, A. Gavdush, A. Alekseeva, P. Nikitin, I. Dolganova, P. Karalkin, A. Khalansky, I. Spektor, M. Skorobogatiy, V. Tuchin, and K. Zaytsev, “Terahertz dielectric spectroscopy and solid immersion microscopy of ex vivo glioma model 101.8: brain tissue heterogeneity,” Biomed. Opt. Express **12**, 5272–5289 (2021). [CrossRef]

**49. **S. Oh, S.-H. Kim, Y. Ji, K. Jeong, Y. Park, J. Yang, D. Park, S. Noh, S.-G. Kang, Y.-M. Huh, J.-H. Son, and J.-S. Suh, “Study of freshly excised brain tissues using terahertz imaging,” Biomed. Opt. Express **5**, 2837–2842 (2014). [CrossRef]

**50. **Y. Zou, J. Li, Y. Cui, P. Tang, L. Du, T. Chen, K. Meng, Q. Liu, H. Feng, J. Zhao, M. Chen, and L.-G. Zhu, “Terahertz spectroscopic diagnosis of myelin deficit brain in mice and rhesus monkey with chemometric techniques,” Sci. Rep. **7**, 5176 (2017). [CrossRef]

**51. **S. Di Meo, J. Bonello, I. Farhat, L. Farrugia, M. Pasian, M. P. Camilleri, S. Suleiman, J. Calleja-Agius, and C. Sammut, “The variability of dielectric permittivity of biological tissues with water content,” J. Electromagn. Waves Appl. (2021). [CrossRef]

**52. **K. Okada, Q. Cassar, H. Murakami, G. MacGrogan, J.-P. Guillet, P. Mounaix, M. Tonouchi, and K. Serita, “Label-free observation of micrometric inhomogeneity of human breast cancer cell density using terahertz near-field microscopy,” Photonics **8**, 151 (2021). [CrossRef]

**53. **I. Dolganova, P. Aleksandrova, P. Nikitin, A. Alekseeva, N. Chernomyrdin, G. Musina, S.-I. Beshplav, I. Reshetov, A. Potapov, V. Kurlov, V. Tuchin, and K. Zaytsev, “Capability of physically reasonable oct-based differentiation between intact brain tissues, human brain gliomas of different who grades, and glioma model 101.8 from rats,” Biomed. Opt. Express **11**, 6780–6798 (2020). [CrossRef]

**54. **A. Ishimaru, *Electromagnetic Wave Propagation, Radiation, and Scattering: From Fundamentals to Applications*, 2nd ed. (IEEE/Wiley, 2017).