Wavefront-modified vector beams for THz cornea spectroscopy

: Terahertz spectroscopy is a promising method to diagnose ocular diseases, where the cornea is typically imaged by Gaussian beams. However, the beam’s mismatch with the cornea’s spherical surface produces a 5-10 % error in analysis. We investigate cornea spectroscopy with wavefront-modified vector beams, reducing the original analysis error to less than 0.5 %. Vector beams are synthesized by our developed 3D Angular Spectrum Method expanded to vector spherical harmonic presentation, allowing wavefront modification and scattering analysis from 100-layer cornea models. We show that wavefront-modified spherical vector beams possess increased accuracy and non-sensitive focusing on cornea spectroscopy compared to the Gaussian beams. Additionally, we investigate wavefront-modified cylindrical vector beams, which show frequency-dependent scattering power arising from s-and p-polarizations. As a result, these beams are unsuitable for cornea spectroscopy, although they have potential for optical force applications. Wavefront-modified vector beams can be applied to spherical target spectroscopy and optical force applications, such as medicine, medical imaging, and optical tweezers.


Introduction
Remote sensing of spherical objects with electromagnetic beams has numerous applications in non-destructive evaluation (NDE) [1], biology [2], and medicine [3].One potential medical application is non-invasive corneal water content sensing for assessing corneal health.Cornea is the transparent front part of the eye covering the iris, pupil, and anterior chamber.The cornea's optical properties are dependent on its water contents, which is ∼ 79 % by mass.The water is distributed as a gradient in the thickness dimension where the top (anterior) is less hydrated than the bottom (posterior surface).The base layers (endothelium and Descemet's membrane) are responsible for maintaining corneal tissue water content (CTWC) by actively pumping out water naturally diffusing in from the aqueous humor.
Diseases, such as Fuch's endothelial dystrophy [4] degrade the water maintenance function leading to increased CTWC.Further, some immunologic responses, such as corneal graft rejection, are typically preceded by, or concomitant with, corneal edema.Later-stage edema is detectable via visual assessment as the hyperhydrated cornea appears "milky" owing to the wavelength-dependent scatter of engorged collagen fibers.However, treatment timed by detection via gross appearance often has limited efficacy.Early and accurate quantification of corneal edema is critical to the management and care of these ocular diseases and detrimental immune responses.
Current diagnostic methods restrict in vivo CTWC quantification to pachymetry, which infers tissue hydration from thickness quantified via ultrasonic or optical backscatter measurements.However, the mapping from the pachymetry-measured central corneal thickness (CCT) to corneal tissue hydration does not account for physiologic variation and returns erroneous, abnormal hydration values for corneal thicknesses off the nominal 580 m [5].In other words, water content is not measured directly but inferred from morphology.
Recent research suggests terahertz (THz) systems are suitable for mapping and directly measuring human cornea water content by quantifying specular THz reflection or coupling to the cornea's lossy longitudinal modes [6][7][8].When referenced to THz band wavelengths ( ∼ 3 -0.3mm), the cornea presents as an optically smooth, wavelength order thickness, spherical shell structure with radially varying dielectric properties arising from the water content gradient in the axial dimension [9].Furthermore, natural physiologic variation in corneal radius of curvature (RoC ∼ 7.8 mm ± 0.5 mm), when referenced to THz band wavelengths, is negligible; thus, the cornea can be considered a perfect sphere with fixed RoC.These properties indicate THz spectroscopy and THz thin film metrology as suitable candidates for direct measurement and quantification of CTWC.
In many reports describing THz illumination of the cornea, the incident field is modeled/assumed as a Gaussian beam in which the RoC of the wavefront matches the spherical surface at the optical axis [10,11].As the focused spot size of the incident field is less than the en face radius of the cornea, the electromagnetic properties of the cornea can be modeled with effective media theories (EMT) [12] to compute the effective permittivity at a given depth, and stratified media theory (SMT) [13] to capture the dielectric gradient and describe the aggregate, frequency-dependent complex reflection coefficient.As the Gaussian beam axis and corneal axis are colinear, and the phase front curvature is matched, the field is considered normal incidence across the entire illuminated area.Optimization routines are then applied to estimate corneal water content and thickness from backscattered THz illumination, assuming the equivalent normal "plane wave incidence on a planar, infinite transverse extent stack." While the planar assumption and subsequent analysis offer reasonable accuracy, the Gaussian beam phase front is hyperbolic and diverges significantly from a sphere within a few confocal distances ( 0 =  2 0 /) from the beam waist  0 , the range over which most corneal sensing is done, see Figure 1.The non-spherical phase front produces spatially varying, non-normal incidence off-apex, resulting in up to 5 -10 % error when analyzing complex reflection with the normal incidence SMT model [10].This paper uses a recently developed 3D Angular Spectrum Method (3D ASM) [14] to synthesize beams that have Gaussian amplitude distributions combined with perfectly spherical phase fronts at the point of corneal incidence.As the cornea shape and dielectric properties are radially symmetric, the spherical coordinate system provides three symmetric orthogonal polarization distributions that define the E-field tangential to the surface and the k vector pointed to the spherical center.The 3D ASM enables the definition of the beam shape coefficients (BSC), and then Mie scattering (cornea size parameter:  = 2/ ≈ 16 -160) is used to compute the broad band complex reflection coefficient for each candidate distribution, and the results are compared to the planar, SMT model.These results determine the optimal polarization orientation to produce a backscatter that most closely resembles the planar model, thus allowing the application of simple plane wave theory to THz spectroscopic data of corneal reflectivity.

Theory
The desired incident fields are synthesized from electric field distributions positioned on the spherical surface (E surf ) using the 3D ASM, presented in detail in [14] and summarized in the Supplemental document.The key concept is to divide the surface area into differential planar elements used as local angular spectrum planes for the 3D ASM.Then, the total field is obtained as a superposition of VSH expanded angular spectrum elements (AS element) to compute the Mie scattering from the layered sphere.This approach allows us to define the ideal electric field and its phasefront on the surface of the cornea, but the field itself can be created elsewhere.
Incident field synthesis from spherical electric field distributions can be obtained simply by defining base vectors for the AS element's local coordinate systems and following the presented procedure in [14].For simplicity, only the base vector presentations and final formulas are presented.

Field synthesis from a spherical surface
Each AS element is presented in a local coordinate system by base vectors (e 1 , e 2 , e 3 ) and an origin at point o(, ).The local beam propagates along the e 3 direction with the electric field polarized in e 1 and magnetic field in e 2 , are obtained from the global spherical unit-vectors (e  , e  , e  ) at a AS element on the spherical surface, see Figure 2. The vector beams are created on the spherical surface by selecting spherical base vectors and the proper propagation axis.The left panel of Figure 3 is termed a spherical vector Beam (SVB) [15,16] with the E field polarized in the  direction and the beam centered at (, ) = (90 • , 0 • ).The middle and right panels are termed cylindrical vector beams (CVB) [17][18][19][20][21] and are centered at (, ) = (0 • , 0 • ).The middle panel CVB is polarized in the radial direction (), and the right panel CVB is polarized in the azimuthal direction ().
(1) Fig. 3. Candidate field polarizations.Note that only the E-field is shown.For all distributions, the E field is tangential, and the k vector is pointed towards the sphere's (cornea's) center The SVB propagates to the negative x-direction and is synthesized from the electric field grid positioned around point  = 90 • .CVBs propagate to a negative z-direction and are synthesized from the electric field grid positioned around point  = 0 • , see Figure 4.The surface electric fields for beam synthesis are Gaussian distributed and defined at a spherical surface as where the arc length from the optical axis is defined as  =  for azimuthal CVB and  =  cos −1 [sin  cos ] for SVB and radial CVB,  is the radius of the sphere and  0 = 1.6 mm is the minimum beam focus.The surface electric field is sampled with Fibonacci distribution to obtain equal AS element spacing and to avoid  = 0 element.Note that these names reflect the similarity of these beams with standard CVBs, but they are not CVBs as there exists no beam waist plane where the polarization is entirely planar.
After following the procedures in article [14], the total electric field presentation is obtained as where r is the location vector for P, and E  (r) is the VSH expanded differential AS element defined as where M 1 emn , M 1 omn , N 1 emn and N 1 omn are the even and odd VSH of the first kind and   emn ,   omn ,   emn and   omn are modified beam shape coefficients (BSCs) including each AS element´s position and orientation information [14].Equations 3 and 4 are explained in more detail in Supplemental documents.Then similar scattered field presentations mapped with the T-matrix method are presented as where E  is replaced with E  , introduced in the Supplemental document.The incident and scattered magnetic fields can be computed similarly by rearranging the BSCs and VSHs, introduced in the Supplemental document as well.

Results
Section 3.1 illustrates the simulation geometry, followed by Section 3.2 introducing SVB, radial, and azimuthal polarized CVB field components on the evaluation plane P1 and propagation plane P2, see Figure 5.Then, Section 3.3 shows the scattering analog between planar SMT and spherical structures with vector beams by 27-layer bandstop filter simulations.Section 4 studies cornea spectroscopy, and Section 4.1 concludes with the effects of vector beam off-focusing.

Simulation arrangement
Analytical formulas for SVBs and CVBs can be derived from Maxwell´s equations using paraxial approximation [15,19,20,22], where numerical aperture (NA) is below NA 0.7 [23].In these cases, the incident fields are defined from the beam waist, preventing wavefront modification at the cornea intersection.In this research, vector beams are synthesized from electric field distributions on the spherical surface using Eq.(3) and Eq. ( 4).This approach creates beams with NA= 0.57 for SVB, NA= 0.77 for azimuthal CVB and NA= 0.89 for radial CVB at the central frequency 300 GHz, where both CVBs are clearly outside of the paraxial approximation range.
Following simulations are realized with the electric field distribution of 45 • total angle on the spherical surface with Gaussian distributed surface field amplitude Eq. ( 2) to maximize the spatial matching between the incident and scattered fields by preventing interference patterns in scattering.Then, incident and scattered fields are computed at plane P1 transverse to the propagation axis at 40 mm from the origin of the sphere and plane P2 along the propagation axis.The evaluation plane P1 location is selected to match the approximated location of the lens in an optical system, where the converging beam can be focused with reasonably small focusing elements, see Figure 5.
As the 3D ASM approximates the beams from the curved surfaces, the error in the incident field synthesis was compared at the evaluation plane P1 to the identical rigorous physical optics simulation at the center frequency 300 GHz of upcoming simulations, with precise matching down to less than a −41 dB average difference in amplitude and less than 1 • average difference in phase.This ensures the appropriateness of the approximation for the study in question.

Spherical and cylindrical vector beams
Wavefront-modified SVBs are a logical candidate as they can be produced by altering the phase distribution of a linearly polarized Gaussian beam.Conversely, CVBs possess a spherically symmetric polarization about the apex and more robust longitudinal modes, which can, theoretically, increase the spectroscopy accuracy of spherical objects.These beams are evaluated on planes P1 and P2 at 300 GHz center frequency in Figures 6-8, showing the Cartesian components of the incident field at the desired optical element location.P1 is a 60 x 60 mm 2 plane at a 40 mm distance from the origin of the sphere.The  2 parameter was computed numerically for all three beam candidates as an initial assessment of the effect that imposing a spherical phase front has on the converging beam.The minimum waist radius and beam divergence angle were computed and used to compute the incident beams quality factor  2 : where  is beams divergence half angle,  0 =1.6 mm is the beam minimum focus radius and  =1 mm is the wavelength.Gaussian beams' 2 is 1.Wavefront-modified vector beams have quality factors The SVB, whose amplitude distribution is TEM 00 Gaussian, yields an  2 1.96; almost twice that of the fundamental mode.The radial and azimuthal CVBs are 3.02 and 2.89, respectively.These results are consistent with observations in [10], namely, optimization on spherical surface offset from the approximate focal point degrades the beam focusing at the focal point.The radial and azimuthal CVBs possess similar total incident field distribution at plane P1.The difference between the beams manifests in the Cartesian components and total field orientation [21].Also, the difference can be seen in the propagation distribution at P2, as the radial polarization includes stronger longitudinal components.All the beams are non-symmetric around the beam waist, as they lack the univocal focal plane.

Scattering analog between planar and spherical structures
To maximize corneal spectroscopic accuracy, the scattering from spherical layered structures should approach the SMT cornea model defined as plane wave scattering from a planar equivalent structure.To explore scattering from these structures and the interaction of the beam with the layers, two different scenarios were considered: (1) concentric spherical shells with lossy cores and (2) concentric spherical shells with lossless cores and, consequently, multiple internal reflections.Additionally, their analogous planar counterparts were considered, see Figure 9.The coupling coefficient K, defined in the Supplemental document, is a selected measure to observe scattering, as it presents scattered energy and spatial matching between the incident and scattered fields.Further, this measure plays a critical role in the cornea´s CTWC and thickness analysis, where the obtained K in the observed frequency range is fitted to match the cornea´s CTWC-and thickness-dependent reflectivity curve from the planar SMT model.
The backscattering behavior from layered spherical structures illuminated by wavefrontmodified SVB and CVBs, along with a an ideal, focused TEM00 Gaussian beam, were studied by reproducing a planar 27-layer bandstop filter to an analog spherical filter structure.Then, the scattering from both filter structures is compared via coupling coefficient .The ideal Gaussian beam was defined such that its on-axis radius of curvature matched the spherical target radius of curvature at a point on the axis more than one confocal distance from the beam waist [24].
The original filter was a 3627 nm thick planar structure with alternating layers of TiO 2 and SiO 2 , with refractive indexes of 2.5 and 1.45, respectively [25].A spherical equivalent structure is scaled to the THz region with a 7.8 mm radius sphere with a lossless and lossy (  = 1 − 0.5) core to prevent internal reflections.Then, the spherical bandpass filter is illuminated in the 200 -300 GHz frequency range by SCBs, CVBs, and by super-confocally focused Gaussian beam presented in [24,26].These results are compared to the plane wave illumination of the original planar filter structures in Figure 10.The simulation shows an increased and well-matched K between the reference planar SMT model and SVB compared to the super-confocal Gaussian beam illumination.This result validates similar scattering behavior between spherical structures illuminated by wavefront-modified SCBs and the SMT model illuminated by the plane wave.In this simulation example, the spherical layers extend the entire radius of the sphere from the surface to the core, allowing the beam to disperse as it propagates to the deepest layers.This phenomenon produces a slight difference between spherical SVB spectroscopy and planar SMT reference.Interestingly, the radial and azimuthal polarized CVBs show decreased matching to the SMT model compared to the SVBs.
In both target cases, the backscatter from the SVB and CVBs produces a closer match to the PW equivalent structure than does the Gaussian beam, thus demonstrating the benefit of phase front matching.The traces in Figure 10 also show the deviation between Wavefront-modified vector beams on spheres, and the planar equivalent is exacerbated at lower frequencies as it is difficult to maintain normal incidence across the interrogated area as the wavelength increases.

Cornea spectroscopy by spherical and cylindrical vector beams
The present study investigates cornea spectroscopy using wavefront-modified SVBs and CVBs, where the corneal central thickness (CCT) was modeled as a 580 -680 m thick, 7.8 mm radius of curvature, spherical shell overlaying a lossy water core.The gradient permittivity of the cornea was modeled with a 100-layer structure using effective medium theory via the Bruggeman model [27].The lossy water core permittivity was calculated with the double-Debye model [28], and the anterior water content (AWC) of the cornea was fixed at 40 %, while the posterior water content (PWC) was varied between 70 % and 90 %.First, the total fields E  = E  + E  of SVB and CVBs from the 100-layer cornea model were simulated in a plane along the propagation axis to illustrate the beam behavior and the equal phase wavefront on the spherical surface, as shown in Figure 11.The coupling K of the cornea model was computed at evaluation plane P1 using four different illumination scenarios within the 200 GHz -400 GHz band: SVBs, radial and azimuthal polarized CVBs, and Gaussian beams with frequency-dependent waists focused at the super-confocal distance [10].The values were compared to the reference values obtained from the SMT cornea model with the same layer thickness and permittivity values, as depicted in Figure 12.The coupling K of the SVBs shows a perfect match with an amplitude less than 0.1 % and a phase less than 0.05 % average difference compared to the SMT model.Gaussian-beam illumination produces a 5-10 % difference for both values.This improved matching in K reduces the fitting errors in the cornea's central thickness and water content (CTWC) analysis.
The coupling of the radial CVB decreased, while the one for azimuthal CVB increased compared to the SMT model.This unexpected phenomenon is confirmed by simulating the scattered power with time-averaged Poynting vectors P  through the plane P1, defined in the Supplemental document.The coupling K and scattered power are also simulated from the homogeneous spheres, which gives the same changes in coupling and scattered power with azimuthal and radial polarized CVBs at lower frequencies, leveling off when moving to higher frequencies.
To understand the changes in coupling and scattering, especially the increased values, the incident fields are simulated at the boundary of the cornea with a frequency of 200 GHz, where the difference in K, with respect to the plane wave (PW), is most significant, see Figure 13.The surface distribution (Eq.2) for synthesizing SVBs and CVBs are initially created with phase-matched electric field AS elements, assuming a unified total electric field on the surface.However, the actual incident field on each element on the surface is a combination of fields from all AS elements, and the polarization of the incident field defines the local interference.As shown in Figure 13, radial CVBs varying phase front and electric field's radial components correspond to the p-polarized oblique angles on the surface.
Similarly, azimuthal CVBs varying phase front and magnetic field's radial components correspond to the s-polarized oblique angles on the surface.SVBs have almost desired tangential electric and magnetic fields on the spherical surface with a small number of radial components corresponding equally to s-and p-polarization, which average polarization effect to zero.These radial field components are frequency-dependent, leveling off at higher frequencies, and the surface fields with 1 THz frequency did not have an affecting amount of radial components.
In conclusion, both CVBs will have varying phase front and longitudinal components on the spherical surface due to the interferences from electric field polarization.These components will affect the oblique angles, creating s-or p-polarized interaction with the surface, reducing or increasing reflectivity.On the other hand, wavefront-modified SVBs maintain the phase front and tangential surface fields on the spherical surface, which enables analog transmission and scattering behavior to SMT models.

Vector beam off-focusing error
Corneal in vivo spectroscopy suffers from significant beam placement and alignment limitation, as the patient eyes are under continuous movement throughout any diagnostic procedure.For this reason, the effect of ±1 mm off-focus is studied by varying the location of E surf at the optical axis.Coupling  and scattered power P  of vector beams and Gaussian beam with focused and off-focused wavefronts are compared in Figure 14, where the cornea is modeled with 580 m CCT 70 % PWC.The results show that SVB maintains scattering behavior through ±1 mm off-focus.The perfect coupling indicates spatial matching between the incident and scattered fields, and evenly distributed s-and p-polarization components do not create frequency-dependent scattered power.However, CVBs are sensitive to beam focusing, which can be explained by the s-and p-polarization components.As the wavefront focus is outside the cornea, the diverging incident beam increases the oblique angles at the surface, leading to increased scattered power.Vice versa, inside-focused wavefronts decrease the oblique angles, leading to fewer changes in scattered power.On the other hand, off-focused Gaussian beams do not have a frequency-dependent scattered power due to the same amount of s and p-polarized components.These components are verified similarly to in Figure 13.Despite the matching scattered power between the Gaussian beam and SVB, the Gaussian beam coupling decreases evenly over the entire frequency range because of the spatial mismatch between the incident and scattered fields.
The results are summarized in Table 1.

Conclusions
Wavefront-modified vector beams (SVBs and CVBs) offer clear advantages in spherical structure spectroscopy compared to nominal Gaussian beams.By matching the wavefront to the corneal surface, the scattering behavior of SVBs approaches stratified medium theory (SMT) reference, which describes the plane wave scattering from an analogous planar structure.This similarity in scattering between SVBs and SMT increases corneal spectroscopy accuracy by closely matching the SMT model, which describes the complex reflectivity according to the corneal thickness and water distribution and the corresponding material parameters over the frequency band.
To demonstrate scattering similarity, the coupling coefficient between the incident and scattered fields from SVBs, radial and azimuthal CVBs, and super-confocally focused Gaussian beams illuminating a spherical 27-layer bandpass filter was computed and compared to SMT reference scattering from a plane wave on an analog planar filter structure.The results showed increased coupling with wavefront-modified vector beams compared to traditional Gaussian beam illumination.
To further validate the proposed method, the gradient permittivity of the cornea was modeled as a 100-layer spherical structure, and the coupling between an incident and scattered fields was evaluated in the 200 GHz to 400 GHz region.The results showed that the SVBs had the closest coupling agreement with the SMT model, with less than 0.5 % amplitude and 0.1 % phase average difference.This indicates that SVBs interact with the multilayered sphere similarly to a plane wave interacting with planar structures, maintaining their beam waist through scattering and producing scattered amplitude and phase that match the SMT model.This improves the accuracy of terahertz frequency corneal thickness measurement and computed tomography using effective medium theory (EMT) and SMT.
Contrary to expectations, the radial CVB showed decreased coupling, while the azimuthal CVB had even higher coupling than the SMT model.This was attributed to the varying scattered power, verified by computing each beam's incident and scattered energy through the evaluation plane.This phenomenon was explained by the orientation of s-and p-polarized electric and magnetic fields on the cornea's surface with CVBs.Radial CVB possesses longitudinal electric field components corresponding to p-polarization and decreased reflectivity.In contrast, azimuthal CVB possesses longitudinal magnetic field components corresponding to s-polarization and increased reflectivity.This phenomenon is frequency-dependent and levels off at higher frequencies, and the surface fields with 1 THz frequency did not have radial components.
Longitudinal components make CVBs more sensitive to focusing errors, as off-focus changes the oblique angles at the spherical surface.Oblique angles increase as the beam´s modified wavefront moves away from the surface, increasing the coupling and scattered power difference.To conclude, CVBs do not match the SMT model, and the error increases with off-focusing, making them impractical in clinical use.
It is important to note that improved coupling to the lossy modes of the cornea in support of THz thin film comes at a cost, as traditional optics cannot produce modified SVBs and CVBs.Future realizations of these beams will require novel optical modalities such as metaoptical lenses [29], metasurface-based optics [30], or specially designed phase plates [31].

Fig. 1 .
Fig.1.The incident beam propagates towards the sphere, and the mismatch between the beam wavefront and the spherical surface is proportional to the distance from the optical axis.

Fig. 4 .
Fig. 4. Electric field grid on a spherical surface, where local electric field, magnetic field, and propagation direction are presented with the (e 1 , e 2 , e 3 ) vectors, marked as red, green, and blue respectively.On a) is CVB base for both polarizations, on b) is SVB base, and on c) is Fibonacci distributed AS elements on a spherical surface.

Fig. 5 .
Fig. 5. Simulation arrangement.a) Incident CVB field and b) incident + scattered CBV fields on both evaluation planes P1 and P2.SVB fields are computed similarly.

Fig. 6 .
Fig.6.Incident SVB with spherically matching wavefront.On the top are Cartesian components at evaluation plane P1 and below the total field on both evaluation planes P1 and P2; see Figure5

Fig. 7 .
Fig.7.Incident radial CVB with spherically matching wavefront.On the top are Cartesian components at evaluation plane P1 and below the radial total field on both evaluation planes P1 and P2; see Figure5

Fig. 8 .
Fig. 8. Incident azimuthal CVB with spherically matching wavefront.On the top are Cartesian components at evaluation plane P1 and below the rotating total field on both evaluation planes P1 and P2; see Figure 5.

Fig. 9 .
Fig. 9. Analog between spherical and planar layered structures.On a) and b) are lossless spherical and their equivalent planar structures.On c) and d) are lossy spherical and equivalent planar structures.

Fig. 10 .
Fig. 10.Coupling coefficients from the spherical layered filter model, where a) presents the lossy core with  = 1 − 0.5 permittivity and b) presents the lossless core scenario.

Fig. 11 .
Fig. 11.Total electric field of SVB and CVBs outside the 100-layer cornea model.a) Cornea model with radially changing permittivity.b) Amplitude (dB), on c) phase (rad) of SVB at xy-plane, d) the amplitude and e) phase of CVB at zy-plane.

Fig. 13 .
Fig. 13.The spherical components and phases of incident fields on the cornea surface.On the left column are electric fields, in the middle are magnetic fields, and on the right are Poynting vector components.The x-axis describes the arc angle from the beam's optical axis.

Fig. 14 .
Fig. 14.On a) are coupling coefficients  with ±1 mm focusing, and on b) are scattered power P  with ±1 mm focusing.

Table 1 .
Averaged  at 200 to 400 GHz band from on and off-focused vector beams and Gaussian beam compared to SMT model.