Metasurface Optics With On-Axis Polarization Control for Terahertz Sensing Applications

Nondestructive testing (NDT) and imaging of materials with unknown properties are important applications of terahertz technology. Different contrast-forming methods are available to obtain diverse information on the imaged region, including intensity, color, phase, and polarization. Polarimetric imaging is relatively under-investigated owing to the difficulty in its implementation with complex optical arrangements. In this article, we present the design and monolithic fabrication of an anisotropic metasurface that incorporates several beam-forming functions in a single layer with high efficiency thus simplifying instrument construction. The probe beam generated by the metasurface consists of an orthogonally polarized pair of collinear Bessel beams that create a well-defined on-axis change of polarization. The metasurface is integrated into a polarimetric imaging microscope for 3-D beam profiling and sensing experiments at 2.52 THz. An analytical model to describe the on-axis polarization variation was found to be in excellent agreement with both finite-difference time-domain (FDTD) simulations and the experimental results. The setup was used to measure the refractive index and diattenuation of homogeneous sample regions in a proof-of-principle application. We anticipate this technology will find future applications in NDT and polarimetric imaging of polymer composites and 3-D profilometry of reflective surfaces.

Whilst many advances have been made in developing suitable technologies to operate in the terahertz band, there is still a considerable need for improvement [1], [13]. Compact high-power sources and sensitive detectors are limiting factors for many applications [3]. Bolometric detection at room temperature is widely used for cameras with large pixel arrays despite the high thermally induced noise floor (E ph (T Hz) ∼ k b T ) [2]. The performance of terahertz imaging systems is further compromised by losses arising from unwanted reflections [14], [15], atmospheric and material absorption [9], [15] and misalignment between individual optical components. Standing waves between constituent components can further deteriorate the performance as a result of unwanted interference [14]. As a consequence of low power sources, lossy optics and low detection sensitivity, terahertz setups often suffer from a low signal-to-noise ratio (SNR) [7].
One promising strategy to overcome the noise limitation of terahertz systems is to develop efficient and highly integrated optical components that enable the construction of compact systems with advanced functions [3], [13]. Metasurface optics can compress several components into a single layer [16], [17], [18], eliminating the problems associated with interference occurring between traditional components. Rigorous optimization with finite-difference time-domain (FDTD) simulations and precise fabrication typically result in high transmission (≥70 %) even for complex optical functions [16], [17], [18], [19]. Metasurface optics may therefore simplify the construction of imaging systems to improve the SNR, unlocking the potential for advanced imaging techniques that would otherwise be unfeasible. It has been previously shown that a metasurface can be designed to analyze the state of polarization (SoP) of an incident terahertz beam [20]. In this work, we explore the capability of metasurfaces to manipulate the SoP into complex 3-D profiles [16], [17], [18], [19], [21], [22] and form a probe beam for sensing applications. The SoP describes the spatio-temporal oscillation of the electric field vector of an electromagnetic wave at a given position. For a wave propagating in z direction, the time-averaged SoP is fully described by the transversal E-field components E x , E y and the phase shift δ between them [23].
A dielectric metasurface typically consists of periodic structures with subwavelength dimensions in one [24] or more directions [25]. The smallest periodic volume is defined as unit cell, which includes the subwavelength structures -often referred to as meta-atoms -that manipulate the properties of light in a well-defined manner. Every meta-atom across a metasurface can be tailored to manipulate transmission properties with subwavelength spatial resolution, including phase and birefringence, to impart complex phase and polarization profiles. The ability to impose two independent phase functions onto an orthogonal pair of polarizations (OPP) has been previously demonstrated [16], [17], [19], [21], [22]. Each OPP consists of two SoP situated at antipodal points on the Poincare sphere. The normalized Stokes parameters (S 0 ≡ 1) form a Cartesian coordinate system on the Poincare sphere, with S 1 = +(−) 1 referring to horizontally (vertically) polarized light, |H (|V ), S 2 = +(−) 1 describing (anti-)diagonally polarized light, |D (|A ), and S 3 = +(−) 1 describing right (left) circularly polarized light, |R (|L ).
These metasurface capabilities can be exploited to generate new types of beam-forming optics based on Bessel beams [19], [21], [26]. Bessel beams exhibit a transversal intensity profile that follows a Bessel function of n th order, which is retained along the optical axis in a non-diffracting manner. By introducing independent phase-encoding onto an OPP, a new class of Bessel beam can be generated that changes its SoP along the optical axis [19], [27], [28], [29]. A probe beam of this type has great potential as a beam-forming optic owing to the beam's ability to retain focus over large axial distances, its self-healing capacity after encountering obstructions [29], [30] and the well-defined correlation between propagation distance z and transversal SoP [28]. A wide range of materials exhibit polarization-dependent properties, such as birefringence or diattenuation, that alter the transmitted SoP profile of such a probe beam. These effects can be used for material identification, characterization, and defect analysis [4], [18].
The terahertz probe beam is composed of a collinear pair of orthogonally polarized Bessel beams to generate an extended focal line with linearly varying SoP along the propagation direction. The monolithically etched metasurface compresses the entire beam-forming function into a single layer. The measured etch depth tolerance of ± 2.3% (±σ) was found to only alter the desired phase shift of each meta-pillar by ±7.1 • to ± 22.8 • depending on the pillar geometry (W x , W y ). The new optical component was integrated into a polarimetric terahertz microscope that was used to demonstrate the retrieval of multiple pieces of sample information, including the refractive index and diattenuation, using a model that accurately described the observed on-axis polarization effects.

A. Design and Operation of the Metasurface
The design of our dielectric metasurfaces is based on rectangular pillars (meta-atoms) with 150 µm pillar height H, whose anisotropic widths W x , W y ∈ [15 µm; 30 µm] are optimized together with the in-plane rotation angle θ ∈ [0 • ; 180 • ] for each lattice point l = (x l , y l ) of a hexagonal lattice with the nearest neighbor distance P of 40 µm, following an established routine for polarization-dependent phase encoding [16], [17], [21], [22], as detailed in supplement section S1. Intrinsic silicon (resistivity >10 kΩ·cm) was chosen for both the substrate and pillars because of its high refractive index of 3.418 refractive index units (RIU) [31] and transparency at the design wavelength λ 0 of 118.8 µm. Fig. 1(a) and (b) shows the fabricated metasurface with its optimized unit cell, respectively. The polar plot of Fig. 1(c) summarizes the simulation results of the unit cell for different pillar geometries. Each data point represents the complex transmission (radial distance | r c | = T x,y , polar angle ϕ c = ϕ x,y ) of a specific pillar geometry (W x , W y ) for incident xand y-polarized light (blue squares and orange crosses, respectively). The imparted phase shifts ϕ x,y span a full phase period of [−π; +π] while the transmission T x,y remained high and uniform at an average value of 86.9% ± 6.4% for both polarizations, close to optimal performance (green circle).
The metasurface was designed to produce an extended focal line in which the SoP changes linearly as a function of the 2-D x-y phase profiles (axicons) generated by the pattern algorithm of supplement section S1 according to (1) to generate different phase profiles ϕ 1,2 with (a) radial phase period g 1 of 0.24 mm for incident |L polarized light and (b) g 2 of 0.4 mm for incident |R polarized light, respectively. The white arrows and dashed circles indicate the radial phase periods g 1,2 and the black dots indicate lattice points l at which the optimized meta-pillars were placed. distance along the optical axis z. This was done by encoding two spatial phase functions ϕ 1,2 (x l , y l ) that are independently excited by an OPP, in our case by incident |L or |R polarized light [27]. All parameters resulting from these two phase functions are subsequently denoted by the subscripts "1" and "2". Both phase functions ϕ 1,2 (x l , y l ) were implemented as axicons with different radial phase periods g 1 < g 2 according to (1), with the radial distance | r l | to the center of the metasurface and x l , y l the Cartesian coordinates of each respective lattice point l. .
(1) Fig. 2(a) and (b) shows the phase shifts (color code) imparted on incident |L and |R polarized light by the selected meta-pillars at each lattice point (black dots) to approximate the phase functions ϕ 1 and ϕ 2 with g 1 and g 2 of 0.24 mm and 0.4 mm, respectively. The pattern algorithm (supplement section S1) accurately approximated both orthogonal phase profiles.
Collective simulations using the software Lumerical FDTD solutions are shown in Fig. 3 for the metasurface of Fig. 2 with radius R MS of 2 mm (9'061 individual pillars). The x-z intensity profiles (y=0) in Fig. 3(a) and (b) confirm that both phase functions ϕ 1 and ϕ 2 work independently to generate different Bessel beams, B 1 or B 2 , upon excitation by incident |L or |R light, respectively. The half-cone angles α 1,2 = sin −1 (λ 0 /g 1,2 ) [27] and non-diffracting focal lengths z max1,2 = R MS / tan α 1,2 of the two Bessel beams B 1,2 are illustrated in Fig. 3(a) and (b) as white dashed lines that accurately describe the simulation results. The black dashed lines describe the transverse diameter of the focal line d f 1,2 = 2.405 · λ 0 /(π sin(α 1,2 )), with d f 1 < d f 2 [30]. The transverse circular area (square area for simulations) around the optical axis spanned by the diameter 2/3 · d f 1 will henceforth be used to define the SoP at the beam center for any given on-axis position z. The SoP along the focal line remained circularly polarized with inverted handedness for both individual beams of Fig. 3(a) and (b).
Upon illuminating the same metasurface with |H polarized light as shown in Fig. 3(c), or any other linear SoP, one can observe a linearly polarized beam center the orientation angle ψ of which changes along the on-axis direction z. In this case, Fig. 3. Collective simulations of a metasurface with radius R MS =2 mm, g 1 =0.24 mm and g 2 =0.4 mm designed to independently generate two different Bessel beams for incident |L and |R polarized light (OPP). The axial x-z intensity profiles at y=0 exhibits two independent Bessel beams, namely (a) B 1 for incident |L and (b) B 2 for incident |R polarized light. (c) Incident |H light excites both Bessel beams, the superposition of which results in a linearly polarized beam center with rotating orientation angle Ψ ∝ z along the propagation direction (black arrows). (d) Inserting a homogeneous slab of dielectric material, e.g., high-density polyethylene (HDPE, between pink dashed lines) into the beam path upon |H excitation alters this superposition as a function of sample thickness t s and refractive index n s (2.0 mm and 1.58 RIU, respectively). Analysis of the (e) orientation angle ψ(z) and (f) ellipticity angle χ(z) at the beam center (data points) provided a detailed analysis of the on-axis SoP with and without the HDPE slab in the beam path. The analytical model of (2)-(4) (lines) is in excellent agreement with the respective simulation results (data points).
both constituent Bessel beams are excited simultaneously. All superimposed electric fields E x1,2 and E y1,2 along the focal line are constant and given by the chosen OPP. However, the on-axis phase shift δ 1↔2 between rays of the two constituent Bessel beams B 1 and B 2 grows linearly with the propagation distance along z according to (2), with the spatial frequencies along the optical axis k z1, 2 This phase shift δ 1↔2 ∝ z translates into a linearly changing SoP along the focal line with a periodicity Λ SoP given by the following equation: The introduction of a homogeneous dielectric slab of sample thickness t s and refractive index n s into the beam path alters the internal half-cone angles β 1,2 = arcsin[sin(α 1,2 )/n s ] according to Snell's law and the spatial frequencies according to k n,z1,2 (n s ) = k 0 · n s · cos(β 1,2 ) for both beams, see Fig. 3(d).
The altered phase shift between B 1 and B 2 at any point on the optical axis z behind the slab is thus given by the following equation: In addition to a standing wave effect caused by Fabry-Perot resonances, Fig. 3(d) clearly shows the decreased halfcone angles β 1,2 within a 2.0 mm thick slab of HDPE (n Lit (HDPE) = 1.58 RIU [15]) between the pink dashed lines. Both Bessel beams B 1 and B 2 were stretched in the axial direction within the sample since β 1,2 < α 1,2 , resulting in an axial displacement of the transmitted rays after the sample region. This effect becomes more apparent when looking at the characteristic angles of the polarization ellipse, particularly the orientation angle Ψ(z) in Fig. 3(e). The observed change of Ψ along the optical axis without slab (red data points) is linear, with a simulated period Λ SoP of 1.38 mm that is in excellent agreement with (3). For the same metasurface but behind the introduced HDPE samples of 0.5 mm (green diamonds) and 2.0 mm thickness (blue triangles), the previously mentioned axial displacement is observed as the shift of the orientation angle ΔΨ ∝ t s that is accurately described by our analytic model using (4). The on-axis ellipticity angle χ(z) of the SoP in Fig. 3(f) remains close to the expected value of 0°until the shorter beam B 1 of Fig. 3(a) becomes increasingly faint beyond ca. 3.5 mm.
The on-axis SoP along one axial period Λ SoP is illustrated on the Poincare sphere in Fig. 4, from which one can infer the role of the metasurface's underlying OPP. For the metasurface design of Fig. 3 that responded to incident |L and |R light (black cubes), the on-axis SoP (red data points) circles around the equatorial S 1 -S 2 plane with an axial frequency of z/Λ SoP . However, a respective metasurface responding to incident |D and |A light follows a circular path along the meridial S 2 -S 3 plane. In conclusion, our metasurface design provides control over the SoP formed along the focal line by choice of its underlying OPP (electric fields E x1,2 and E y1,2 ), whereas the oscillation length Λ SoP can be chosen via the difference in radial phase periods Δg = g 1 − g 2 (opposite rotational sense for g 1 > g 2 and g 1 < g 2 ). The presented metasurface design enables various polarimetric applications in terahertz sensing and imaging by simply scanning the metasurface along the z-axis with respect to the object under study, as will be shown in the experimental results section.

B. Metasurface Fabrication and Evaluation
We developed a single-step photolithographic process to dry etch the metasurface layer into high-resistivity silicon substrates (25 mm x 25 mm, 535 µm thick) using an ICP-DRIE Bosch method based on C 4 F 8 /SF 6 plasma chemistry [32]. A detailed description is given in supplement section S2. Fig. 5(a) shows a scanning electron micrograph of the pillar's sidewall profile exhibiting a close to vertical angle of 88.5 • -89.5 • . Fig. 5(b) shows a plan view of a larger area of the metasurface. Optical interference profilometry (Bruker, ContourGT-X) was used to obtain a 3-D profile of the pillars and their surrounding etch floor, as shown in Fig. 5(c). A variation was observed in the etch depth owing to an effect known as aspect-ratio dependent etching [32]. The measured 3-D data was analyzed as a histogram over the depth dimension. The average etch depthH = (z pillar − z f loor ) is determined reliably with its standard deviation σ H = σ 2 pillar + σ 2 f loor from a Gaussian fit to the histogram using the two peak positions z pillar and z f loor and their respective variances σ 2 pillar and σ 2 f loor . The measured Simulations of the pillars were rerun with adapted pillar height ΔH and substrate thickness of ± 5 µm (≈ 1.5 σ H ) and ± 10 µm (≈ 3 σ H ) to gauge the impact of the etch depth deviation on the device performance, as detailed in supplement section S3. All pillar geometries exhibited a linear relation (R ≥ 0.95) between the resulting phase and height differences Δϕ x,y ∝ ΔH [25] with a slope ranging from 2.1 • /µm to 6.7 • /µm, or on average (4.1 ± 1.0) • /µm. Fig. 5(d) shows the central X cross-sections through the encoded axicon phase profiles ϕ L , ϕ R along the white dashed line shown in Fig. 5(c). For both ϕ L and ϕ R , the phase shifts encoded by the individual pillars (black data points, H = 150 µm) accurately form the desired conical phase fronts (black lines). Taking a worst-case etch depth deviation of ± 5 µm -equal to 1.5 σ H -into account resulted in the red error bands for the imposed phase profiles.

C. Experimental Setup for Polarimetric Imaging Microscopy
The experimental setup to use the beam-forming metasurface for polarimetric imaging of the beam profile and NDT applications is shown in Fig. 6. An optically pumped continuous wave laser (FIRL 295, Edinburgh Instruments) delivered a collimated Gaussian beam (1/e radius ca. 5.5 mm) of up to 120 mW at 2.52 THz, which was |H polarized. A non-polarizing beam splitter with power meter in its side arm was introduced to account for fluctuations in the output power between individual images. A pair of wire grid polarizers allowed for dynamic attenuation within the measurement arm by adjusting the angle of WGP 1 while leaving WGP 2 fixed to transmit |H . A custombuild telecentric objective that projected the 7.7x magnified detection plane z det on to the camera's sensor (see supplement section S4) was placed behind the sample. A rotatable quarter wave plate (QWP), a fixed analyzer in |V transmission and a bandpass filter were integrated in the objective's infinity space. The combination of the QWP and fixed polarizer allows for the measurement of the SoP [23]. The bandpass filter is used to increase the SNR at the terahertz camera (INO MicroXCam-384-THz).  Fig. 7 shows experimental results obtained using the fabricated metasurface of Fig. 5 with pattern radius R MS of 11 mm and radial phase periods g L and g R of 0.72 mm and 0.39 mm that generate |L and |R polarised Bessel beams, respectively. The half-cone angles α L and α R were predicted to be 9.5°and 17.8°, creating a theoretical oscillation period of the on-axis SoP Λ theo of 3.45 mm using (3).

III. EXPERIMENTAL RESULTS
The transverse polarization profile can be generated pixel-bypixel from eight individual images with incrementally rotated QWP by 22.5 • using the setup of Fig. 6, see supplement section S5 for a detailed description. Such polarization-resolved images of the probe beam are shown in Fig. 7(a) and (b), depicting the characteristic angles Ψ(x, y) and χ(x, y) of each pixel's polarization ellipse, respectively. The transverse profiles of both 0 th order Bessel beams form equidistant concentric rings around the beam center in Fig. 7(a) and (b), with the white and black dotted circles representing d f,L and d f,R , respectively. The beam center exhibited uniform Ψ and χ close to 0 • , corresponding to |H light. The circular mean of the pixel values within the beam center (diameter of 2/3 · d f,R ) was calculated as reliable and reproducible SoP for each respective scan position z i along the focal line. The residual |L and |R components in Fig. 7(b) are a result of the steeper half-cone angle, creating a more tightly focused beam center d f,R < d f,L for the |R wavefront. The more tightly focused beam center for |R caused a slightly imbalanced superposition (I R > I L ) of the two constituent Bessel beams B R and B L with some residual ellipticity at the beam center (I res = I R − I L ) and the radially alternating rings of predominantly |L or |R components seen in Fig. 7(b). The residual center ellipticity χ can be exploited to assess the quality of polarimetric images since it provides information on potential misalignments between the constituent beams B R and B L , and whether the maximum scan range z max (B R ) of the shorter beam has been surpassed (I R I L ). In contrast, the SoP's orientation angle Ψ is exclusively defined by the superposition of equal circular components (I L ≡ I R ) and their phase difference δ L↔R . Ψ is thus robust against perturbations caused by intensity variations in I res or misalignment and is therefore well suited to be a measurement channel in the following sensing experiments. Fig. 7(c) shows the progression of the SoP's orientation angle Ψ along the focal line obtained in a calibration scan (blue data points). The observed linear course is in excellent quantitative agreement with our model of (3), (4) (dashed magenta line) with Λ exp of 3.57 mm very close to its theoretical value of 3.45 mm. The setup's capability to analyze a sample's diattenuation was tested using a randomly oriented wire grid polarizer. The axial SoP beating observed in the calibration measurement changed into a sine-shaped on-axis intensity beating with Ψ(z) ≈ const. (orange diamonds) after the polarizer was introduced into the beam. The known orientation of the polarizer's transmission axis of Ψ T = 25 • (green line) was measured to be 26.2 • ± 3.7 • by averaging all scan positions (orange dashed line).
Finally, Fig. 7(d) demonstrates the ability to measure the refractive index n HDP E of a sample of the known thickness (t HDP E = 5.65 mm). The on-axis orientation angle Ψ(z) at the beam center exhibits a shift ΔΨ = (Ψ Cal (z) − Ψ HDP E (z)) averaged to 78.3 • ± 6.0 • (violet data points) after transmission through the HDPE sample compared to the calibration, as expected from the simulations shown in Fig. 3(e) and the analytical model of (2)-(4). This SoP shift ΔΨ(z) ≈ const. was used together with the known sample thickness t HDP E to retrieve the measured refractive index (cyan data points) using the analytical model. The measured refractive index of 1.55 RIU ± 0.06 RIU was obtained by averaging over all data points of this graph and is close to the literature value of 1.58 RIU (red line) [15]. Note that the measurement sensitivity S RI = ΔΨ/(t s · [n s − 1]) (the SoP change in response to change in optical path length) is set by the metasurface design parameters. Increasing Δg = g 1 − g 2 will increase α 1,2 and β 1,2 (n s ) hence δ 1↔2 (z, t s , n s ) and thus the measured signal ΔΨ. The dynamic measurement range Δn max = [π /(S RI · t s )] − 1 will decrease as a tradeoff, limiting the range within which refractive indices exhibit an unambiguous shift ΔΨ ≤ π. The measurement sensitivity S RI of the metasurface employed in Fig. 7

IV. CONCLUSION
We have presented an innovative concept for a terahertz beam forming metasurface that exploits the SoP as an information channel. The versatile metasurface design is capable of controlling the transmitted phase-and polarization-front using anisotropic dielectric pillars to form a probe beam. The component was monolithically fabricated by etching into high resistivity (>10 kΩ·cm) silicon to yield an experimental transmission efficiency of up to 76.0 ± 0.4%. Detailed evaluation of the fabrication tolerances demonstrated a practically insignificant degradation of the metasurface's optical performance. The result is a single compact optic that minimizes the refraction and diffraction problems associated with constructing a multicomponent system [21], [27], [28], [29].
An analytical model was developed to predict the axial progression of the SoP of the resulting probe beam from its constituent pair of orthogonally polarized Bessel beams. The model was shown to be in good agreement with detailed simulations and experimental data.
The metasurface was embedded as a beam-forming condenser lens in a polarimetric imaging microscope and both the probe beam's transverse profile and the on-axis change of the SoP were found to closely match theoretical predictions. The microscope was used in sensing applications including the measurement of the diattenuation and refractive index of a homogeneous sample. The resulting instrument for refractive index sensing is the first of its kind to operate at terahertz wavelengths and showed improved measurement sensitivity and dynamic range when compared to previous work at a wavelength of 532 nm [21]. Compared to established techniques such as terahertz timedomain spectroscopy (THz-TDS), our approach can provide rich polarimetric information on the sample, e.g., diattenuation and birefringence. Respective polarimetric measurements using THz-TDS typically require two measurements with a rotated sample by 90°, thereby only analyzing two orthogonal linear SoPs rather than the complete polarisation ellipse. Furthermore, our approach allows for the use of large pixel arrays to image a wide field-of view (FOV), whereas THz-TDS imaging typically requires slow lateral pixel-by-pixel scanning of the sample. The image acquisition time of the presented setup is limited by the slow SoP measurement (rotating QWP method [23]) that can be accelerated to as little as 1.2 s when combined with a previously published metasurface polarimeter [20]. In terms of data processing, we expect our technique to be easily automated and compare favorably with THz-TDS with respect to data volume and processing when considering images of similar FOV since 8 intensity measurements per pixel are processed for our technique compared to the temporal information of each pixel in THz-TDS. THz-TDS has the advantage of providing spectral sample information. Although our demonstrator setup used a monochromatic light source, we expect from previous studies that the type of metasurface presented here should maintain good operation over a bandwidth of more than 1 THz [20]. Consequently, the presented metasurface design could also be used as an optic in a THz-TDS system with the potential to add further functionality to existing instruments.
Future work should explore further optimization and extension of both beam forming and setup capabilities to enable polarimetric imaging of heterogeneous samples. Furthermore, a new imaging instrument is envisioned that would be capable of analyzing the reflection from a sample with the incident probe beam. Such a setup could retrieve a 3-D surface profile from polarimetric information and the diattenuating or birefringent behavior of the sample simultaneously. Such an instrument could in the future have a role in NDT of polymers and resin-based composites.
Yash D. Shah received the B.E. (Hons.) degree in