Open Access
Add to BibSonomy
Abstract
Coupling to longitudinal modes of thin spherical shells, under Gaussian-beam illumination, was explored with a theoretical method based on Fourier-optics analysis and vector spherical harmonics and was scrutinized with an experimental setup. For the theory part, the illumination frequency band was fixed between 100–600 GHz and the outer spherical shell radius of curvature and thickness are 7.5 mm and 0.5 mm, respectively. The shell material was either the lossless cornea or an aqueous effective media representing the cornea. Six different beam-target strategies were introduced being potential candidates for maximum coupling. Two dispersion-tuned beam ensembles with strongly frequency-dependent phase center location have been created with a fixed incident beam 1/e radius and radius of curvature called forward strategies. These computations of different alignments were continued with four beam ensembles of frequency-invariant phase center, constructed from fits to experimental data, oriented at four different axial locations with respect to the spherical shell center of curvature, they are called reverse strategies. Coupling efficiency for all strategies was calculated for different targets including perfect electrical conductor (PEC) sphere, PEC core covered by a cornea loss-free layer and cornea. All scattering strategies contrasted to scattering from equivalent planar targets as a reference with maximum coupling. The results show that, under an ideal calibration, forward strategies are a closer approximation to the plane-wave condition for the cornea. An experimental setup was assembled to explore the simulation approach in a frequency range between 220 GHz to 330 GHz. Two different quartz samples with permittivity of 4.1 were mounted on a water core, acting for a cornea. The first and second quartz radius and thickness were 7.5 mm and 0.5 mm and 8 mm and 1 mm, respectively. An adequate agreement between theory and experiment was confirmed. A particle optimisation swarm algorithm was applied to extract the thickness and permittivity of quartz from the measured back-scattered field for reverse strategies.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The sub-millimeter wave and THz frequency sensing of the cornea leverages the layered tissue structure for corneal water content and thickness quantification. Changes in corneal tissue water content (CTWC) and corneal central thickness (CCT) are correlated with human eye diseases and disorders. Existing clinical measurement approaches are not sufficiently accurate for the early detection of these conditions. A promising method for non-invasive, in situ tissue characterization, is THz spectroscopy. The cornea is $~ 0.6$ mm thick and is bounded by air and an optically thick body of water on its anterior (outer surface ) and posterior (bottom inner) surface, respectively. At THz frequencies, this structure presents as a lossy thin-film lying atop a lossy half-space. Resolution of the cornea’s lossy longitudinal modes via frequency-domain reflectometry in a band sufficiently low (e.g. $220$ GHz - $330$ GHz) for significant penetration allows simultaneous estimates of CTWC and CCT. However, since the cornea is spherical, an efficient coupling to longitudinal modes nominally requires normal incidence across the interrogated area and thus a converging spherical phase front whose curvature matches the corneal surface curvature [1–3]. However, as the average corneal radius of curvature (RoC) is $\sim 7.5$ mm and the midband free-space wavelength at $100$ GHz - $600$ GHz is $\sim 1$ mm, approximate phase-front matching occurs near the beam confocal point where the distinctly non-spherical phase-front curvature is rapidly changing and leads to low longitudinal modes coupling.
A typical, non-contact measurement system block diagram is shown in Fig. 1. A transceiver optical co-locates a source and detector to enable normal incidence illumination and reception. The output beam is collimated and focused onto the cornea and back-scattered reflection is referenced to some calibration standard reflection to compute broadband, complex sample reflectivity. Fits between the measured data and a model are used to estimate water content or water content gradients. Additionally, spectra can be used to estimate thickness and axially varying hydration simultaneously.
Fig. 1. Scheme of the problem. Modelling cornea as a planar layered structure under plane wave illumination as well as a spherical layered structure under Gaussian illumination is shown.
The nominal output beam profile of a THz system is TEM$_{00}$ Gaussian and different strategies regarding the placement of focused beam waist with respect to corneal RoC have been evaluated [4]. Restated from above, while the Gaussian beam phase front radius of curvature can match the cornea RoC on-axis, there is no configuration that achieves perfect matching across the field of view. Nevertheless, spectra are almost exclusively analyzed with Fresnel’s equations at normal incidence informed by effective media permittivities to model depth-dependent water content and stratified media theory to model the depth-dependent and thus provide an aggregate reflection coefficient of the stack. Implicit in this analysis is a planar thin film on semi-infinite planar half-space, with plane wave illumination. The radius of curvature mismatch at the cornea’s air-anterior interface and posterior-aqueous humour interface and resulting walk-off losses within the cornea itself cannot be accommodated.
In summary, THz spectroscopy of cornea is performed via focused Gaussian beam illumination on lossy hemispherical shell enclosing a lossy homogenous aqueous core while analysis assumes planar illumination on a planar, infinite transverse extent lossy thin film on top of a semi-infinite aqueous half space. The above-listed constraints raise two interrelated questions about obtaining maximum coupling: ($1$) where should the corneal center of curvature (CoC) be positioned relative to beam geometry? ($2$) Is there a frequency dependence on Gaussian beam parameters that further helps to maximize coupling? More succinctly, given a set of constraints on the beam, how does one maximize coupling to the cornea’s sub-millimeter wave longitudinal modes? It’s worthy to say maximum coupling is considered as convergence to plane-wave on planar surface coupling due to the perfect match of plane wave’s phase-front with the planar surface.
To best address, the above questions, a method based on Fourier optics (FO) and vector spherical harmonics (VSH) is used to computationally explore longitudinal coupling in structures resembling cornea. The VSH is the vector solution of the wave equation in spherical coordinates. In 1908, Gustave Mie was one of the pioneers who used VHS to address the incident and scattered field from a sphere illuminated by a plane wave [5]. Later, Davis modified the theory for Gaussian beam illumination [6]. The early proposed method was quite intricate and cumbersome. In 1993, Khaled et al. [7] presented a method using Fourier analysis [8] to model Gaussian beam as the angular spectrum of plane waves, and for computing scattering coefficients, they employed the T-matrix method [9].
In fact, the T-matrix relates the scattering coefficients of a target illuminated by a plane wave with the scattering coefficients of the target illuminated by a Gaussian beam. To compute the scattering coefficients of layered sphere under the plane wave illumination many approaches were presented [10,11,12,13], here Yang [13] algorithm was employed. This classical method will serve as a new application in our paper allowing us to scrutinize the cornea as a multilayered spherical structure. FO allows for a closed-form expression of the steady-state scattering solution (multiple reflections within the corneal layer) while VSH provides a convenient representation of the fields with respect to the target geometry. The goals of the paper are to address different beam target strategies which are potential candidates for maximum coupling and investigate these strategies theoretically and experimentally for spherical targets. Additionally explored was the efficiency of using a metal sphere as a calibrator [14] in an experimental set-up.
2. Background and Gaussian-beam analysis
At this point, the constraints on the illuminated beam parameters are explored, addressing the questions in the introduction, whether there is a frequency dependence on Gaussian-beam parameters assisting coupling enhancement. Also, to determine which beam-cornea alignment are the most likely optimal strategies to explore.
2.1 Cornea as a spherical scatterer and spherical cavity
Interestingly, a cornea can be considered as a scatterer in a wide frequency range. A plot of the size parameter range for cornea in the THz band is shown in Fig. 2(a) with wavelength along the horizontal axis and particle radius along the vertical axis. The family of oblique lines is parameterized by a fixed size parameter, $ka=2\pi n_0/\lambda$, where $a$, $k$, $\lambda /n$ are respectively, the radius, the wavenumber, and the wavelength in the medium in which the particle resides ($n_0 = 1$ for air refractive index). The fixed size parameter contours of $0.0002$, $0.2$, and $2000$ serve as approximate thresholds between the varying scattering regimes. The Mie scattering is a proper candidate to describe the scattering behavior of particles with size parameters lying in the range $0.2< ka < 2000$.
The vacuum wavelength range $3$ mm - $0.3$ mm corresponding to the $100$ GHz - $1$ THz band is indicated by the gray shaded vertical box. The nominal range for human cornea RoC, $7$ mm - $8$ mm is denoted by the thin purple horizontal shaded rectangle. Its intersection reports the approximate corneal size parameter range subtended by the anterior segment. The size parameter of the posterior segment was estimated by subtracting the nominal CCT range from the anterior segment RoC and computing the wavelength in the complex aqueous corneal medium. Bruggeman’s effective media theory [15] is used to compute the aggregate permittivity of a mixture of $60\%$ collagen and $40\%$ free water.
Fig. 2. (a) Cornea scattering region, (b) stability analysis of the cornea as a spherical cavity, coupling efficiency computed by ABCD matrix for beam illumination (c) on the apex, and (d) on the center of the cornea.
These intracorneal size parameters are labeled "anterior segment" and "posterior segment" in Fig. 2(a). The dispersion arising from the liquid water component is evident in the curvature of the shaded region. At $100$ GHz, both the anterior and posterior segments lie comfortably in the Mie scattering regime. At $1$ THz, the posterior segment is quite close to the approximate geometric scattering limit and, depending on the illumination profile, indicating maybe GO is sufficient to approximate the expected back-scatter.
Moreover, the cornea can be treated as a spherical cavity. It implies certain constraints on the beam. Only particular ranges of the cavity outer radius $R_1$, inner radius $R_2$, and the distance between them $L$, produce stable resonators. An unstable cavity will increase the beam size without limitation, consequently, it will get larger than the cavity size and will be lost completely.
A stability analysis of the cornea for nominal RoCs and thicknesses are displayed in Fig. 2(b) where $g_1 = 1-L/R_1=0.9333$ and $g_2 = 1-L/R_2=0.9286$. The unstable region $g_1g_2>1$ is indicated by the shaded region and contours of constant $g_1g_2$ by the dotted lines. For clarity, both RoC sign pairs are plotted; ($+R_c, -R_c+t$) and ($-R_c, +R_c-t$). In a region about the corneal apex, the anterior and posterior surfaces are concentric and thus $g_1g_2 = 1$ for all combinations of $R_c$ and $t_c$. The corneal thickness increases towards the periphery which can be described as an anterior segment RoC increase with respect to the posterior segment RoC and thus a divergence from concentricity. This combination is unstable and is represented by the overlap between the $g_1g_2$ curves and shaded unstable region.
Thus, the corneal cavity is susceptible to beam walk-off which reduces the interference between the primary reflection at the air/cornea interface and multiple reflections from paths through the cornea. Two examples of beam walk-off are shown in Fig. 2(c) using ABCD matrix computations [16]. The cornea was modeled via the effective media theory described above but only the real part of the permittivity was applied. Ray transfer matrices were used to compute the parameters of the beam scattered from the anterior segment and then the beams scattered after one, two, and three round trips through the cornea. The plots show the coupling coefficients (Eq. (11), [4]) as a function of frequency between the primary reflection beam and the round trip beams. The superior coupling of the RoC-matched beams compared to the incident waist radius beams suggests that the phase front matching will enhance interference and thus aid in the coupling and extraction of longitudinal modes.
2.2 Gaussian-beam analysis
The problem geometry definitions are displayed in Fig. 3. A converging beam is an incident, from the left, on the cornea. The incident plane is labeled plane $1$ and the incident beam transverse radius and RoC are denoted $\omega _1$ and $R_1$, respectively. The beam waist is labeled $\omega _0$ and located at plane $0$ where the parameters of plane $0$ and plane $1$ are related via free-space propagation in the absence of the cornea. The physical distance between plane $0$ and plane $1$ is labeled $d_{01}$ and the distance between the corneal CoC and plane $0$ is $z_0$.
Two beam ensemble definitions were considered where the incident beam RoC was fixed such that $R_1= R_c \forall$ frequencies. The first was termed "forward" (shown in Fig. 3) where $\omega _1$ was set at some constant value and then ray transfer matrices and complex beam parameters were used to compute $d_{01}$ and $\omega _0$ as a function of frequency:
Fig. 3. Geometry of the problem, cornea modelled as a layered sphere under Gaussian beam illumination by two different beam ensemble definitions including forward and reverse strategies.
The second approach was termed "reverse" or "varied confocal distance" (shown in Fig. 3) where $\omega _0$ as a function of wavelength is defined and the standard Gaussian beam equation for axially dependent RoC was used to find $d_{01}$ and then $\omega _1$. Note that $d_{01}$ is the solution to a second order polynomial thus two solutions are possible. These are labeled $d_{01,NF}$ and $d_{01,FF}$ in Eq. (3) and (4), respectively where the $d_{01,FF} > d_{01,NF}$ for $z_c<R_c/2$. The iF subscript stands for $i=N,F$ indicating FF far-field and NF near-field, respectively:
The $d_{01}$ and $\omega _0$ in the forward direction over a range of $\omega _1$ are reported in Fig. 4(a,b) for $100$ GHz - $1$ THz. In both plots, the ($f$, $\omega _1$) pairs that produce waist radii that violate the paraxial limit $\omega _0< \lambda /2$ are eliminated. Further (f, $\omega _1$) pairs that produce $d_{01} < Z_c$ are denoted. Input beam radii $\omega _1 = 2.1$ mm and $3.1$ mm are denoted by the white dotted lines and are considered in the next section.
Fig. 4. (a) Forward $d_{01}$ and (b) forward $\omega _0$ over a range of $\omega _1$ is plotted. (c) Reverse $d_{01,NF}$ and (d) reverse $\omega _{1,NF}$ over a range of $\omega _0$ is plotted.
Parameters spaces for $d_{01,NF}$ and $\omega _{1,NF}$ computed in the reverse direction are displayed in Fig. 4(c,d) with the ($f$, $\omega _0$) pairs in violation of the indicated paraxial approximation. The large area in the upper right-hand corner of the plot correspond to ($f$, $\omega _0$) pairs where $Z_c>R_c/2$ and thus complex valued $d_{01}$ has been set aside. Inspection of Eq. (4) shows that the contour lines in the $d_{01,NF}$ plot of Fig. 4 correspond to contours of constant $z_c$. In other words, if one defines $\omega _0$ such that $Z_c$ is invariant to frequency then, as evidenced by Eq. (3), $d_{01,NF}$ and $d_{01,FF}$ are invariant to frequency and thus the beam RoC evolves equally along the axis for all spectral components.
2.3 Gaussian-beam strategies
According to previous discussions, the final candidate strategies are displayed in Fig. 5 where the left column shows the evolution of the beam radius and beam waist location with respect to the corneal geometry and the middle column shows the evolution of the beam RoC superimposed on the corneal anterior and poster surface locations. The right column reports the spectral dependence of beam waist $\omega _0$, incident beam radius $\omega _1$, and beam waist location $Z_0$.
Fig. 5. Definition and orientation of Gaussian beams and their relationship to corneal geometry for (a-c) $S_1$, (d-f) $S_2$, (g-i) $S_3$, (j-l) $S_4$, (m-o) $S_5$ and (p-r) $S_6$ strategies.
Strategy $1$ ($S_1$) fixes the input beam radius to $\omega _1=2.1$ mm yielding a diameter similar to that of modern ultrasound pachymeter probes [17]. Additionally, $S_1$ allows analysis of illumination where the low-frequency incidence occurs in the near-field region of the beam and the high-frequency incidence in the far-field. This is evident by the black trace denoting the confocal point location (Fig. 5) and requires significant dispersion visible by the beam waist locations reported in the left and right columns.
The beams in strategy $2$ ($S_2$) were also defined in the forward direction with $\omega _1 = 3.1$ mm fixed for all frequency. This ensemble locates the incidence location in the far-field region for all beams but is still sufficiently small to avoid $\omega _0< \lambda /2$. The $S_2$ beam waists are slightly smaller than $S_1$ but the dispersion (variation in $Z_0$) is significantly reduced.
The reverse analysis was utilized for strategies $3$ - $6$ ($S_3$ - $S_6$) with the confocal distance fixed at $z_c = 2.62$ mm. This confocal distance is below the $R_c/2$ threshold and yields physically realizable, paraxial approximation compatible, providing a set of beams for the frequency range $100$ GHz - $1$ THz.
The $S_4$ places the phase front match in the super confocal (far-field) region of the beam as evidence by the beam overlapping beam RoC plots in Fig. 5. The beam waists are collocated and the waist radii are all larger than $\lambda /2$ although they approach the limit at $100$ GHz. The $S_3$ is the same beam ensemble as $S_4$ but places the phase front match in the subconfocal (near field) region of the beam. The beam radius on the cornea is significantly smaller ($\omega _1 \sim \omega _0$) which should reduce phase front mismatch error in the beam periphery but the beam RoC magnitude is rapidly increasing for increasing $z$ instead of decreasing suggesting a substantial mismatch at the posterior surface.
Strategies $5$ and $6$ ($S_5$, $S_6$) were evaluated for comparison to two strategies commonly reported in the literature. The $S_5$ places the beam waists at the corneal CoC. This is the typical arrangement for imaging via a Gaussian beam telescope optical train and produces phase fronts that are slightly larger in RoC than the corneal RoC. The $S_6$ places the beam waist at the corneal apex thus mimicking the common approach reported by many groups using THz time-domain spectroscopy. The beam radius on target is minimized at the cost of a significant RoC mismatch. These strategies are labeled "reference" throughout this work and are iterations of the reverse strategies.
3. Angular spectrum based analysis
An approach based on Angular Spectrum, for calculating the coupling efficiency of sub-millimeter wave illumination on the cornea, was presented in [18] which utilized the methodology described in [7]. This method is applicable for any Gaussian beam incident on a spherical surface and is valid for spheres with RoC in order, or larger than the illumination wavelength. The advantage of this method compared to PO and full-wave approaches are the ability to solve directly the steady state scattering of both homogeneous and coated spheres. This ability significantly reduces computation time as many issues such as "ray splitting" at dielectric interfaces and finely discretized surface current densities can be avoided. Thus, it is possible to explore several illumination strategies on the scattering profile of the cornea when it is modeled as a spherical shell encasing a lossy dielectric sphere.
The applied approach first decomposes the Gaussian beam to its plane-wave spectrum representation and then expresses these components as VSH [7]. The incident Gaussian beam is represented as:
For the calculation of coefficients $a^c_n$ and $b^c_n$, Khaled in [7] and [19] utilized the algorithm introduced by Toon and Ackerman [10]. In this work, the algorithm employed by Yang [12,13] was applied to the coefficients calculations. The equations for the algorithm is addressed in Appendix B. Electromagnetic analysis for a layered sphere with more than one layer is the advantage of the Yang algorithm.
Figure 6 shows the field distribution of incident Gaussian beam and back-scattered field from cornea computed with 5 and 6, in a plane perpendicular to the optical axis. The plane dimension is $20$ mm by $20$ mm. The source beam was polarized in the $x$-$z$ direction. The cornea was modeled as a single-layered spherical shell encapsulating a homogeneous, pure-water sphere. The water permittivity was obtained by the double-Debye model [20], and the shell (cornea) modeled with the effective medium theory via the Bruggeman mixing model [15]. The corneal shell consisted of $60\ \%$ water and $40\ \%$ collagen with dispersion-free and real permittivity $2.9$. Core radius was $7.5$ mm and shell thickness was $0.5$ mm.
Fig. 6. Field distribution of source and back-scattered field of the cornea is simulated with Eqs. (5) and (6) at 220 GHz. Source amplitude and phase (rad) in $x$, $y$, and $z$ direction is shown in (a,d), (b,e), and (c,f), respectively. Amplitude and phase (rad) of the back-scattered field from the cornea in $x$, $y$, and $z$ direction is shown in (g,j), (h,k), and (i,l), respectively. Field distribution is calculated on a plane perpendicular to the optical axis with dimension $20$ mm $\times$ $20$ mm. The source beam is assumed to be polarized in the $x$ and $z$ direction.
Fig. 7. PEC sphere coupling efficiency magnitude comparison is presented for (a) Forward strategies, $S_1$, and $S_2$, (b) reverse strategies, $S_3$, and $S_4$, and (c) reference strategies, $S_5$, and $S_6$. Plane wave reflection from equivalent PEC planar is indicated by black-dotted line and radius of sphere equals $R_c=7.5$ mm.
Source fields distribution plots based on Eqs. (5) and (6) are almost consistent with the analytical field plot of the Gaussian beam. Since the source beam is polarized in $x$ and $y$ directions, the $E_y$ field amplitude (Fig. 6(b)) is in micro order, much smaller than the $E_x$ and $E_z$ amplitude (Fig. 6(a,c)). The back-scattered field amplitude diminishes about $~ 40\ \%$ in $x$ and $Z$ directions. Likewise, if the cornea is considered as a planar structure illuminated by a plane wave, the reflection amplitude would be also around $40\ \%$ at $220$ GHz. The back-scattered field phase is shifted $180^{\circ }$ as expected.
4. Physical-optics analysis
The six different strategies were simulated with an in-house developed physical-optics (PO) script to verify the results of the presented Fourier-optics method. The two-way propagation from an emitting screen to the homogenized PEC sphere and back was simulated with the Gaussian-beam parameters introduced in section 2. In PO, the field outside the radiating aperture was calculated as in [21], [22]:
The PO simulated electric field was oriented to ensure no shadowing between the surfaces occurred with an edge taper sufficient to limit the spill-over loss. The coupling coefficient from the PO simulation is overlaid with that from the Fourier-optics method, validate it for homogeneous targets such as PEC sphere. For all strategies, the coupling coefficient with FO method and PO differ by less than $1\%$ across the $220$ GHz to $330$ GHz. The coupling coefficient from the physical-optics simulation is consistently less than that from the Fourier method and is likely due to the remaining spill-over loss necessary to avoid shadowing in this geometry.
5. Results and discussion
The coupling efficiency of the six different strategies addressed in section 3 was compared for different targets and compared with the plane wave reflection from a multilayered planar surface as the maximum coupling efficiency case. It was computed with Fresnel equations and superimposed on the Gaussian-beam equations to evaluate divergence from the plane-wave condition. The cornea was modeled as described in section 3 and Fig. 6.
To compare different strategies for various objects, coupling efficiency was defined as Eq. (11). It calculated the coupling between the back-reflected field and the incident field in a reference plane located at $z=-40$ mm. This distance is chosen according to experimental setup. The integration range for $x$ and $y$ starts from $-4\omega _1$ to $4\omega _1$ and $^*$ denotes the complex conjugate. For accurate calculations of coupling efficiency, the field at the edge of the reference plane should be vanishingly small, thus the reference plane was increased until this criteria was met.
where $\mathbf {E}^{i}$ and $\mathbf {E}^{s}$ are computed by Eq. (5) and Eq. (6). The coupling efficiency is defined as the absolute value of coupling coefficient $|CE|$ and phase of $CE$ is defined as $\arctan (\frac {image(CE)}{real(CE)})\frac {180^{\circ }}{\pi }$.Normally, the number of modes ($N_{stop}$) used in CE calculations were determined by the introduced equation in appendix B. To achieve a more robust result, each simulation accuracy was checked in terms of convergence of coefficients to a specific value. For each target, $N_{stop}$ were reported in a table. Simulation time increases exponentially by increasing the number of modes. The frequency band $100$-$600$ GHz was chosen regarding the simulation time and stability of the algorithms as well as constraints introduced in section 2.
5.1 Homogeneous PEC sphere coupling efficiency magnitude
In this subsection, the coupling efficiency of the previously described designs in section 2. for the PEC sphere was calculated using the proposed methodology in section 3. The Gaussian-beam illuminating PEC was assumed to scatter more energy back at the screen than any other material of equal RoC and thus served as a reference/calibration target for the further exploration of PEC coated with the lossless cornea, and cornea simulations.To compute the scattering coefficients of a PEC sphere illuminated by a Gaussian beam Eq. (8) was used. For applying the T-matrix method it was first necessary to compute the scattering coefficients of the PEC sphere when illuminated by a plane wave. Here, the approach in [23] and presented in appendix C was applied. The PEC sphere radius was set to $R_c=7.5$ mm and the beam waist radii locations were distributed to match strategies $S_1$ to $S_6$ as described in section 2. The number of modes to reach enough accuracy for the PEC sphere is reported in Table 1. The simulation time for each strategy was about $17$ hours on UPC computational cluster.
Table 1. Number of modes for PEC sphere and coated PEC sphere simulations
According to Fig. 7, the forward strategies $S_1$ and $S_2$, with dispersion set to achieve frequency wavefront matching, behave similarly and increase from $97.84 \%$ to $99.92 \%$ and from $97.83 \%$ to $99.89 \%$. Both reverse strategies $S_3$ and $S_4$ coupling efficiencies behave almost in the same way and range from $96.46 \%$ to $100 \%$ and $96.77 \%$ to $99.94 \%$ across the band. The reference strategies $S_5$ and $S_6$ behave likewise and unlike the other strategies decrease over the frequency band from $95.71 \%$ to $94.59 \%$ and $96.08 \%$ to $94.59 \%$, respectively.
Overall, for homogeneous PEC sphere being a target, $S_1$ comes closest to the plane-wave condition although it behaves closely to $S_2$. Forward strategies reveal higher coupling compared to reverse strategies ($\sim 1.1 \%$) and reverse strategies display more coupling compared to reference strategies ($\sim 0.4-5.5 \%$ across frequency band).
5.2 Lossless cornea sitting on a PEC core coupling efficiency
Lossless dielectric spherical shells backed by PEC were simulated to explore the effect of spectral phase-front variation and mismatch. A target was constructed of $0.5$-mm thick lossless cornea spherical shell encapsulating a PEC sphere of $7$ mm RoC. The layer refractive index was set as the real part permittivity of the cornea computed by effective media theory [15]. The coupling efficiency was computed with Eq. (11). Then, the calibrated coupling efficiency was computed by normalizing the PEC-backed lossless cornea shell coupling efficiency with the coupling efficiency from a $7.5$ mm PEC sphere (Fig. 8). The scattering coefficients of a coated PEC were addressed in appendix C. The number of modes considered to reach stability for coated PEC was similar to PEC sphere modes in Table 1 and simulation time was the same.
Fig. 8. Calibrated lossless cornea sitting on a PEC core coupling efficiency magnitude and phase comparison are presented for (a,d) Forward strategies, $S_1$, and $S_2$, (b,e) reverse strategies, $S_3$, and $S_4$, and (c,f) reference strategies, $S_5$, and $S_6$. Plane wave reflection from equivalent layered planar is indicated by black-dotted line. The shell permittivity is set according to real part of cornea permittivity and its thickness is $0.5$ mm sitting on $7$ mm RoC PEC sphere.
This operation mimics experimental calibration routines [4] and enables analysis of the absolute phase following deconvolution of the free space path influence on complex coupling angle. The PEC sphere calibrated lossless cornea shell coupling efficiency magnitude and phase are displayed in Fig. 8 where magnitude and phase of forward strategies $S_1$ and $S_2$, reverse strategies $S_4$ and $S_5$, and reference strategies $S_5$ and $S_6$ are plotted in panels (a), (b) and (c), respectively. All trends are referenced to the equivalent plane-wave condition.
Interestingly, as illustrated in Fig. 8, with calibrating all strategies, phase behaviors were almost consistent with plane wave condition, implying a proper phase matching. In fact, they matched perfectly if plane wave phase figure shifts about $10$ GHz toward low frequencies. These results suggest that good results can be obtained with any most of the strategies provided the target reflection is calibrated with reflections from an appropriately sized, metallic sphere. Moreover, from the magnitude aspect, forward strategies provided better coupling between the incident and back-scattered fields than reverse ones, however, $S_4$ featured closer to $S_1$ and $S_2$ (less than $1 \%$ difference). It could be said $S_3$ and $S_5$ were also close enough to plane wave condition (less than $\sim 2.2 \%$ difference) and certainly, $S_6$ was far from proper coupling ($1.5-8.8 \%$ deviation).
The simulation results suggest that the lossless cornea shell is acting as a lens and improving the wavefront match to the $7$ mm RoC inner PEC surface relative to the reference $7.5$ mm RoC PEC sphere. Non-sequential ray-tracing simulations of a Gaussian beam in Zemax OpticStudio also suggest that the presence of a lossless dielectric layer on top of the PEC (hemisphere) sphere might slightly improve the coupling. In Fig. 9, the power distribution on the detector is reported for $S_1$ at two different frequencies, $200$ GHz and $600$ GHz. At $200$ GHz the root mean square (RMS) spot radius is $3.40$ mm for the coated sphere and $3.33$ mm for the PEC sphere. At $600$ GHz the RMS spot radius variation is much more constrained: $3.13$ mm for the coated sphere and $3.12$ mm for the PEC sphere. This trend agrees with Fig. 8(a) where the coupling shows a peak at $200$ GHz.
Fig. 9. Incoherent irradiance for $S_1$ strategy in $200$ GHz for (a) PEC sphere and (c) coated PEC sphere and in $600$ GHz for (b) PEC sphere and (d) coated PEC sphere, indicating higher RMS of coated sample compare to calibration sphere, resulting above $1$ values of coupling efficiency amplitude in Fig. 8. Panel (e) supports the idea by a histogram of normalized frequency of rays in each angle, emphasizing on loss-free shell assists the incoming beam to come closer to the optical axis (lower angles).
In a quasi-optical, mono-static ($S_{11}$) measurement, the maximum signal is obtained when the beam is retro-directive, i.e. the scattered beam mirrors the illumination beam. In the case of a spherical target, this can be visualized by discretizing the wavefront as a spatial collection of converging rays. The maximum signal is achieved when all converging rays are normal to the spherical surface. However, in this case, the converging beam is Gaussian with a non-spherical phase front thus normal incidence is achieved only on-axis. The loss-free shell assists the incoming beam to come closer to the optical axis (as shown in Fig. 9(e)), that’s why the above one values for coupling efficiency appeared in Fig. 8.
5.3 Calibrated cornea coupling efficiency
The approach described in [18], and outlined in Appendix B, yielded the magnitude and phase coupling between the incidence and back-scattered electric fields for cornea structure. The number of modes used for simulation is reported in Table 2 and simulation time was $25$ hours on UPC computational cluster. As mentioned earlier, the best strategy for performing an experiment that matches with the theory developed in previous works [4] is the one that Gaussian-beam illumination on a spherical surface acts such as a plane-wave illumination on a plane resembling a maximum coupling. Magnitude and phase of relative calibrated coupling efficiency for cornea is plotted in the Fig. 10. In this context, relative calibrated means subtract calibrated cornea coupling efficiency from equivalent plane wave condition coupling efficiency.
Calibration with PEC sphere mirrors current experimental methodology where the reflectivity of a phantom or the ex-vivo cornea is normalized by a spherical reference reflector of approximately equal RoC. The reflectivity is obtained on an absolute scale when the optical path is deconvolved from the reflection by the normalization. These results support the previous hypothesis in the literature that, under the right alignment conditions, Gaussian-beam illumination on cornea approximates plane-wave illumination on planar stratified media.
Fig. 10. Relative calibrated cornea coupling efficiency magnitude and phase comparison are shown for (a,d) Forward strategies, $S_1$, and $S_2$, (b,e) reverse strategies, $S_3$, and $S_4$, and (c,f) reference strategies, $S_5$, and $S_6$.
Table 2. Number of modes for cornea simulations
Surprisingly, these strategies behave almost the same as a plane-wave illumination on the planar surface and show a maximum of $2.3 \%$ magnitude deviation in low frequency for strategy $S_6$. Similarly, for calibrated phase coupling efficiency (indicated in Fig. 10(d-f)) of preceding strategies attributes almost the same as plane-wave illumination on the plane, with limited deviation at low frequency for strategy $S_3$ (a maximum of $1.2^{\circ }$ deviation).
These results suggest that, under ideal alignment and calibration conditions, the efficiency of coupling to corneal longitudinal modes is nearly indistinguishable amongst the six different strategies. This behavior is almost similar to the loss-free shell results, although the lossless cornea shell deviates from plane wave condition more than cornea case, indicating that lossy cornea steady state reflectivity minimizes walk-off losses and any illumination condition, regardless of wavefront match/mismatch is sufficient for proper coupling as far as equal surface geometry is feasible for calibration.
6. Experimental exploration
In order to explore the theory was introduced in section 3. and was used for investigation in previous sections, a setup (as shown in Fig. 11(a)) is implemented. Gaussian beam telescope introduced in [4] and shown in Fig. 11(b) is applied to illuminate quartz dome sitting on a water core (resembling a cornea). It is shown previously that a back-scattered measured field from a $0.5$ mm quartz sitting on an air core acts as a plane wave illumination on plane [24]. Here, two sizes of quartz are considered in this experiment, quartz size $1$, $Q_1$ with diameter $15$ mm and thickness $0.5$ mm, and quartz size $2$, $Q_2$ with diameter $16$ mm and thickness $1$ mm. These dimensions are chosen according to real cornea size and their permittivity is around the real part permittivity of the human cornea (around $4.1$). The quartz dome is set on a holder which is mounted on mechanical arms which move in x, y, and z-direction, illustrated in 11(a). The back-scattered field is measured in a frequency range from $220$ GHz to $330$ GHz. As shown in Fig. 11(c), a metal ring is used to hold the quartz dome on the water core and an absorber ring (Fig. 11(d)) covers it to prevent unwanted scattered fields from metal. Figure 11(e) presents the confocal distance of the output Gaussian beam in the quasioptical system. The output beam waist radius of the quasioptical system is measured and fitted with the fundamental Gaussian-beam model. The fit waist radii in E- and H-planes are used to calculate the confocal distance:$z_{c,E}=\pi \omega _{0,H}^2/\lambda$, $z_{c,H}=\pi \omega _{0,H}^2/\lambda$. Due to the slightly different E- and H-plane beam width of the feed horn and polarization effects in the quasioptics, the output beam is astigmatic and the confocal distance is determined separately for the two planes.
Fig. 11. (a) Measurement setup, (b) Gaussian beam telescope illuminate sample, (c) holding quartz dome on a water half-space by a metal ring, (d) use ring absorber to eliminate unwanted back-scattering from metal holder, (e) confocal distance of the output Gaussian beam, (f) OCT image of $Q_1$, and (g) $Q_2$.
Alignment was initiated with a course $2$D transverse (x-y) sweep of the sample about the approximate focal axis of the Gaussian beam telescope. The location corresponding to peak collected back-scattered intensity was identified and then a finer $2$D scan was performed about the peak to further refine the peak location. Next, the sample was fixed on the peak location, and the target was swept 20 mm in the vertical (z-axis) direction over a distance of $20$ mm. Complex $S_11$ were acquired in steps of $0.1$ mm for a total of $201$ axial points. For data analysis, the reflection from the sample at each point was Fourier transformed to the time/space domain, time-gated and transformed back to the frequency domain. This procedure performed repeated for PEC spheres of radius $8$ mm and radius $7.5$ mm and for the water backed quartz dome phantoms.
Examples of range the Fourier transformed data are shown in the Fig. 12(a) and 12(b) for the $7.5$ mm RoC PEC sphere and $7.5$ mm, $0.5$ mm thick water backed dome respectively. The vertical axes corresponds to the approximate physical distance of the target with respect to a reference plane located in the beam path and shows the axial range traversed by the z-axis during scanning. The horizontal axes show the calculated range from the extender output flange to the target as computed with the Fourier transform of the complex $S_{11}$ data. In Fig. 12(a) the (r,z) location of the peak corresponds to the subconfocal point as reported in [4] and the superconfocal point and waist location can be found by traversing the (r,z) space by the known separation distances. The same analysis was applied to the phantom data in Fig. 12(b). Once the locations corresponding to strategies $3$, $4$, $5$, and $6$ are identified in these data, the quartz domes back-scatter can be calibrated with the PEC sphere data. An example of this is shown in Fig. 12(c).
Fig. 12. Examples of range of the Fourier transformed data for (a) the $7.5$ mm RoC PEC sphere, and (b) $7.5$ mm, $0.5$ mm thick water backed dome and, (c) Normalized $S_{11}(t)$ as function of $r$. Purple line shows the spherical mirror calibration data from where the strongest return came (arguably close to the subconfocal distance indicated by green line), sample location which results in best fit to the theoretical is shown with blue line and its pick is in yellow.
In Fig. 13, the measured back-scattered field of the sample is depicted (dashed lines) and compared with VSH simulation results (solid lines). All the plots are calibrated by the same size PEC sphere, both for simulation and measurement. Figure 13(a) and (b) displays the amplitude and phase of quartz size $1$ (Diameter $15$ mm and thickness $0.5$ mm) and Fig. 13(c) and (d) demonstrates the amplitude and phase of quartz size $2$ (Diameter $16$ mm and thickness $1$ mm). There is an adequate agreement between simulation and measurement in different reverse illumination strategies, representing realistic strategies $S_3$, $S_4$, $S_5$, and $S_6$), in order to illuminate a sample as $S_1$ and $S_2$ a custom lens with engineered spectral dispersion is required. Resonances and peak locations in the measured spectrum are consistent with the theory. Also, as the Gaussian beam focus moves from the center of the sample ($|Z_0|=0$ mm) to the apex of the sample ($|Z_0|=7.5$ mm), a slight shift toward the higher frequency is observed both in the measurement and simulation. The shift is more detectable in thick quartz and may be due to the Gouy phase shift and corresponding superluminal velocity which have not been corrected for in this analysis. This result confirms that, with correct calibration, all six illumination strategies produce similar results. Note that the thicker, fabricated quartz dome was $0.96$ mm thick and the thinner one was $0.5$ mm. Both had an estimated $\epsilon \sim 4.0$. These values were used in the VSH simulation. The extraction of quartz dome permittivity and thickness were performed with particle swarm optimization (PSO). The PSO algorithm searches for the best fit to a plane wave model by minimizing a merit function. Strategies $S_3$, $S_4$, $S_5$, and $S_6$ were investigated representing realistic strategies. The merit function combined two terms: the phase error and amplitude error. The used merit function for each scenario is:
Fig. 13. Measurement data are compared with VSH simulation for reverse strategies for (a,b) quartz size $1$ amplitude and phase, and (c,d) quartz size $2$ amplitude and phase. As the beam focuses closer to apex a slight shift toward higher frequency is detected. The shift is more evident in thicker quartz. This is in agreement with simulation results.
Fig. 14. (a) Extracted thickness and (b) extracted permittivity for $Q_1$ are obtained by PSO fitting method considering strategies $S_3$, $S_4$, $S_5$, and $S_6$. (c) Extracted thickness and (d) extracted permittivity for $Q_2$ are obtained by PSO fitting method considering strategies $S_3$, $S_4$, $S_5$, and $S_6$.
7. Conclusion
A theory based on Fourier-optics analysis and VSH was employed to explore the coupling efficiency between the incident and back-scattered fields from homogeneous and coated spherical targets at different beam-target alignments. In this method, Gaussian beam approximated by an ensemble of plane waves and scattering coefficients were calculated via the T-matrix method.
PEC spheres are targets for calibration in different simulations. PEC coupling efficiency magnitude in six different strategies was explored. The $S_2$ alignment shows the highest matching with plane-wave illumination on a planar surface. Also, the influence of $+0.2$-mm RoC change was investigated. The $S_6$ strategy shows the highest sensitivity (close to $3.3\%$) to RoC discrepancies.
To explore the effect of spectral phase front variation and mismatch, coated PEC structure was studied. Due to the lossless layer in this assembly, walk-off loss takes part in the coupling, and beam profile influences the calibrated coupling efficiency magnitude, although phase behavior remains insensitive to the illumination profile and only shows $10$ GHz shift with respect to low frequencies than the plane wave condition. The $S_2$ design revealed the highest consistency with plane wave on plane illumination.
Surprisingly, the calibrated coupling efficiency magnitude for $S_1$ - $S_4$ strategies (in coated PEC) is above one. The reason is the refraction of the Gaussian beam while passing through the loss-free layer, as an incident beam propagates closer to the optical axis, the beam reflection enhances for certain kinds of Gaussian beams. This is in striking contrast to the stratified medium theory and theories considering plane waves - the specifics of the Gaussian beam and target RoC need to be considered carefully before applying the stratified medium model.
Angular spectra and vector spherical harmonic analyses were used to investigate the corneal coupling efficiency. The cornea was modeled as a layered media on top of the aqueous core and illuminated by a Gaussian beam. It was shown that engineering the sub-millimeter wave beam such that incident beam $1/e$ radius fixed to $\omega _1=3.1$ mm and radius of the Gaussian beam matches with $7.5$ mm cornea radius at each frequency (strategy $S_2$) leads to the highest coupling among other strategies.
Illumination profiles defines the coupling efficiency behavior unless the corneal coupling efficiency is calibrated the coupling efficiency of PEC sphere of equal size. In all strategies, coupling efficiency magnitude and phase behave according to the plane-wave illumination on the planar half-space. It implies that the cornea is lossy enough that walk-off loss does not influence coupling efficiency.
An experimental setup is assembled to validate the VSH theory. The experimental back-scattered spectrum is compatible with theory, peak and resonance location are the same for different size of domes and a shift toward higher frequency observed both in experiment and simulation. A PSO analysis is applied to extract the thickness and permittivity of quartz from measured back-scattered field for different realistic strategies.
Appendix A. Gaussian beam coefficients derived by Fourier analysis
For a polarized Gaussian beam in $x$ and $z$ direction, after a cumbersome calculation, the incident beam coefficients in Eq. (5) are obtained as:
In the above equations, the $\Pi =\Pi _{mn}=\frac {P^m_n(\cos \theta )}{\sin \theta }$ and $\tau =\tau _{mn}=\frac {d}{d\theta }P^m_n(\cos \theta )$ are auxiliary functions which can be obtained by recursion relation [9]:
The first two starting values of $\Pi _{mn}$ are generated by the closed-form expressions when $m\neq 0$:
On the other hand, the $A_{\theta _{xy}}=TG\sin \theta _x\sin \theta _y$ where $T=e^{-ik(\textbf {s}.\textbf {v})}=e^{-i(k_xx_0+k_yy_0+k_zz_0)}$ indicates the relocation of waist radius location relative to center of sphere and $G=e^{-(\frac {ks\omega _0}{2})^2}=e^{-\frac {\omega ^2_0}{4}(k^2_x+k^2_y)}$. Besides, $PW=e^{ik(\textbf {s}.\textbf {r})}=e^{i(k_xx+k_yy+k_zz)}$ illustrates a plane wave. The vectors $\textbf {s}$, $\textbf {v}$, and $\textbf {r}$ are $\textbf {s}=\cos \theta _x\mathbf {i_x}+\cos \theta _y\mathbf {i_y}+\cos \theta \mathbf {i_z}$ , $\textbf {v}=x_0\mathbf {i_x}+y_0\mathbf {i_y}+z_0\mathbf {i_z}$, and $\textbf {r}=x\mathbf {i_x}+y\mathbf {i_y}+z\mathbf {i_z}$. The wave vector defines as $\textbf {k}=k_x\mathbf {i_x}+k_y\mathbf {i_y}+k_z\mathbf {i_z}$ and $k=|\textbf {k}|=[k_x^2+k_y^2+k_z^2]^{1/2}=2\pi /\lambda$ is the wave number in the wavelength $\lambda$. For any arbitrary $\textbf {k}$, $\theta _x$, $\theta _y$, $\theta$, and $\phi$ are defined the way that the components of $\textbf {k}$ are:
Appendix B. Coated sphere illuminated by a plane wave
The scattering coefficients of a coated sphere illuminated by a plane wave is presented by [13] and summarized in [12]. The method starts with solving the Helmholtz equation in spherical coordinate and fulfillment of the boundary condition at each layer. The scattering coefficients are derived as below:
Next equations determine the logarithmic derivatives of the Riccati–Bessel functions, $D^1_n(z)$ and $D^3_n(z)$, the ratio $Q^l_n$, $\psi _n$ and $\zeta _n$. A thorough explanation can find in [13] and [12].
As mentioned in section 3, to compute the scattering coefficients of a coated sphere while illuminated by a plane wave, a method different from [7] and [19] was utilized. Khaled et al. used Toon and Ackerman [10] algorithm while in here the algorithm described by Yang [13] is applied. The advantage of the Yang method is providing the possibility of analysis of a coated sphere with more than one shell.
Appendix C. PEC and coated PEC sphere illuminated by a plane wave
The scattering coefficients of a sphere illuminated by a plane wave are computed by solving the vector-wave equation in spherical coordinates and evoking the boundary condition (Mie theory [11]). In the case of a PEC sphere, the electric field at the PEC boundary was forced to zero and the simplified form of scattering coefficients for a PEC sphere are obtained as:
With the mutual approach, solving a vector-wave equation, applying boundary condition and forcing the electric field to be zero inside the PEC core, the simplified scattering coefficients for a coated PEC sphere are obtained as:
The PEC core radius $r_1$ is coated with a layer with outer radius of $r_2$ and refractive index of $m_2$. The Riccati-Bessel functions are $\chi _n(z)=-zy_n(z)$, $\psi _n(z)=zj_n(z)$, and $\zeta _n(z)=zh_n^{(1)}(z)$ in which $j_n(z)$, $y_n(z)$, and $h_n^{(1)}(z)$ are Bessel function of the first, Bessel function of the second kind and Hankel function of the first kind, respectively. The primes denote the differentiation with respect $z$.
Funding
Agencia Estatal de Investigación (PID2019-107885GB-C31/AEI/10.13039, PRE2018-084326); Academy of Finland (327640).
Acknowledgments
We would like to thank Prof. Ari Sihvola for his valuable advice to understand and implement the theory of problem.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
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.
References
1. Z. D. Taylor, J. Garritano, S. Sung, N. Bajwa, D. B. Bennett, B. Nowroozi, P. Tewari, J. W. Sayre, J. Hubschman, S. X. Deng, E. R. Brown, and W. S. Grundfest, “Thz and mm-wave sensing of corneal tissue water content: In vivo sensing and imaging results,” IEEE Trans. Terahertz Sci Technol 5(2), 184–196 (2015). [CrossRef]
2. Z. D. Taylor, J. Garritano, S. Sung, N. Bajwa, D. B. Bennett, B. Nowroozi, P. Tewari, J. Sayre, J. P. Hubschman, S. Deng, E. R. Brown, and W. S. Grundfest, “THz and mm-wave sensing of corneal tissue water content: electromagnetic modeling and analysis,” IEEE Trans. Terahertz Sci Technol. 5(2), 170–183 (2015). [CrossRef]
3. S. Sung, S. Selvin, N. Bajwa, S. Chantra, B. Nowroozi, J. Garritano, J. Goell, A. D. Li, S. X. Deng, E. R. Brown, W. S. Grundfest, and Z. D. Taylor, “Thz imaging system for in vivo human cornea,” IEEE Trans. Terahertz Sci Technol. 8(1), 27–37 (2018). [CrossRef]
4. A. Tamminen, S. V. Pälli, J. Ala-Laurinaho, M. Salkola, A. V. Räisänen, and Z. D. Taylor, “Quasioptical system for corneal sensing at 220-330 ghz: Design, evaluation, and ex vivo cornea parameter extraction,” IEEE Trans. Terahertz Sci Technol 11(2), 135–149 (2021). [CrossRef]
5. G. Mie, “Beiträge zur optik trüber medien, speziell kolloidaler metallösungen,” Ann. Phys. 330(3), 377–445 (1908). [CrossRef]
6. L. W. Davis, “Theory of electromagnetic beams,” Phys. Rev. A 19(3), 1177–1179 (1979). [CrossRef]
7. E. E. M. Khaled, S. C. Hill, and P. W. Barber, “Scattered and internal intensity of a sphere illuminated with a Gaussian beam,” IEEE Trans. Antennas Propag. 41(3), 295–303 (1993). [CrossRef]
8. J. W. Goodman, Introduction to Fourier Optics (Roberts, 2005).
9. P. W. Barber and S. C. Hill, Light Scattering by Particles: Computational Methods (World Scientific, 1990).
10. O. B. Toon and T. P. Ackerman, “Algorithms for the calculation of scattering by stratified spheres,” Appl. Opt. 20(20), 3657–3660 (1981). [CrossRef]
11. C. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley Science Paperback Series, 1998).
12. U. P. O. Pe na, “Scattering of electromagnetic radiation by a multilayered sphere,” IEEE Trans. Antennas Propag. 180, 2348 (2009). [CrossRef]
13. W. Yang, “Improved recursive algorithm for light scattering by a multilayered sphere,” Appl. Opt. 42(9), 1710–1720 (2003). [CrossRef]
14. F. Zarrinkhat, M. Baggio, J. Lamberg, A. Tamminen, I. Nefedova, J. Ala-Laurinaho, E. E. M. Khaled, J. M. Rius, J. Romeu, and Z. Taylor, “Calibration alignment sensitivity in corneal terahertz imaging,” Sensors 22(9), 3237 (2022). [CrossRef]
15. A. Sihvola, Electromagnetic Mixing Formulas and Applications (Institute of Electrical Engineers, 1999).
16. P. F. Goldsmith, Quasioptical Systems: Gaussian Beam Quasioptical Propagation and Applications (IEEE Press, 1998).
17. B. Lackner, G. Schmidinger, S. Pieh, M. A. Funovics, and C. Skorpik, “Repeatability and reproducibility of central corneal thickness measurement with pentacam, orbscan, and ultrasound,” Optom. Vis. Sci. 82(10), 892–899 (2005). [CrossRef]
18. F. Zarrinkhat, J. Lamberg, A. Tamminen, M. Baggio, J. Ala-Laurinaho, E. E. M. Khaled, J. M. Rius, J. R. Robert, and Z. Taylor, “Fourier analysis of submillimeter-wave scattering from the human cornea,” in EuCAP (2021), pp. 1–5.
19. E. E. M. Khaled, S. C. Hill, and P. W. Barber, “Light scattering by a coated sphere illuminated with a gaussian beam,” Appl. Opt. 33(15), 3308–3314 (1994). [CrossRef]
20. J. Xu, K. Plaxco, S. Allen, J. Bjarnason, and E. Brown, “0.15-3.72 thz absorption of aqueous salts and saline solutions,” Appl. Phys. Lett. 90(3), 031908 (2007). [CrossRef]
21. J. Kong, Electromagnetic Wave Theory (John Wiley and Sons, 1986).
22. S. J. Orfanidis, Electromagnetic Waves and Antennas (Rutgers University New Brunswick, 2002).
23. R. Ruppin, “Scattering of electromagnetic radiation by a perfect electromagnetic conductor sphere,” J. Electromagn. Waves Appl. 20(12), 1569–1576 (2006). [CrossRef]
24. F. Zarrinkhat, J. Lamberg, M. Baggio, A. Tamminen, J. Ala-Laurinaho, E. E. M. Khaled, J. M. Rius, J. R. Robert, and Z. Taylor, “Experimental exploration of longitudinal modes in spherical shells at 220 ghz – 330 ghz: applications to corneal sensing,” in CLEO/Europe-EQEC, (2021), p. 1.
