Complex vectorial optics through gradient index lens cascades

Graded index (GRIN) lenses are commonly used for compact imaging systems. It is not widely appreciated that the ion-exchange process that creates the rotationally symmetric GRIN lens index profile also causes a symmetric birefringence variation. This property is usually considered a nuisance, such that manufacturing processes are optimized to keep it to a minimum. Here, rather than avoiding this birefringence, we understand and harness it by using GRIN lenses in cascade with other optical components to enable extra functionality in commonplace GRIN lens systems. We show how birefringence in the GRIN cascades can generate vector vortex beams and foci, and how it can be used advantageously to improve axial resolution. Through using the birefringence for analysis, we show that the GRIN cascades form the basis of a new single-shot Müller matrix polarimeter with potential for endoscopic label-free cancer diagnostics. The versatility of these cascades opens up new technological directions.

The diattenuation value can be readily obtained from the second to fourth elements in the first row of a MM, that is the elements m12, m13 and m14 respectively, as shown in Eq. (5). The retardance R is reconstructed from the trace of , Eq. (6), while the orientation of optic axis ranging from − to radians is calculated according to Eq. (7). The depolarization properties are included in the in the bottom right 33 matrix ∆ of the matrix , which is shown in Eq. (8). , and in Eq. (8) are the eigenvalues of ∆ , and is a matrix composed of the eigenvectors of ∆ . ∆ and ∆ in Eq. (8) are linear/circular depolarization values.
= cos ( ) − 1 , = tan , We found that 1) the measured experimental MM as well as the MMPD parameters matched well with the simulation counterparts; 2) the MMPD results reveal that the birefringence structure of the GRIN lens behaves equivalently to a spatiallyvariant waveplate array with a) linear retardance that gradually increases along the radial direction and b) fast axis directions that vary azimuthally (gradually changing from -π/2 to 0 then to π/2 radians as the azimuthal angle changes from 0 to π radians, then a similar variation when the azimuthal angle changes from π to 2π radians). The fast/slow axis direction distributions are shown in Fig. 1b in the main article as well.
In Supplementary Figures 2d to 2g, we further show four different SHWP based GRIN lens cascades including in Supplementary Figure 2d -the one used in the main article (as Fig. 3b and 3c). Besides the close correspondence between simulation and experimental results, we also found that the MMs of these GRIN lens cascades have significant off-diagonal asymmetry which is of interest in MM research, as well as related sample information analysis 11 . Specifically, we extracted the circular anisotropy coefficient 12 of these four MMs, shown in Supplementary Figure 2h, in which we use grey arrows to indicate the gradient of the circular anisotropy 13 . This implies the system can be used to demonstrate spin-Hall effect (SHE) of light 13 . Both of these phenomena deserve further exploration.
We propose here that the GRIN lens cascade can also find use as a special device for other research areas, such as the derivation of quantitative MM parameters from exotic materials (e.g. chiral characteristics 14 ), biomedical research (multilayer birefringence and/or diattenuation characterization 11 ), or studies on complex spin-orbital interaction (SOI) processes of light 13 . These all take advantage of the unique property of the GRIN lens -its equivalence to a spatially variant waveplate array, which can encompass all combinations of retardance and fast axis direction.  . The blue double-sided arrow is slow axis direction (extraordinary axis), and the orthometric red counterpart is the fast axis direction (ordinary axis) of that arbitrary point. The black and grey focus points at distal side of the GRIN lens represent focus splitting between the e ray (white solid line) and the o ray (white dotted line), which is introduced by both the birefringence property and the sinusoid ray trace. (b) (i) the section of a GRIN lens shows the refractive index and the birefringence profiles, (ii) shows two curves of these two parameters along the radius direction from centre point to the boundary of GRIN lens, (iii) shows an incident light on point A( ′, ′, ′), on the front surface of GRIN lens, referenced by z axis. (iv) shows the lateral profile of one arbitrary section of GRIN lens, which includes a ray trace passing through point B. is the interior angle between the ray and the extraordinary axis (blue double-sided arrow in (a)) at point B, which is the complementary angle to , which is the interior angle between the ray and the direction of the z axis.
Supplementary Figure 3a illustrates the refractive index (violet coloured gradient), birefringence (yellow coloured gradient) as well as one arbitrary ray trace (white solid line) when light propagates through the GRIN lens. Supplementary Figure 3b (i) shows the refractive index profile and birefringence profile of a GRIN lens section, and (ii) shows sketches of these two parameters along the radius direction from the centre to the boundary (refer to the lens from Femto Technology Co. Ltd., G-B161157-S1484). All the calculations below are under the meridional plane approximation 16 .
For a general ray within the GRIN lens, the distance from the GRIN lens axis is denoted by , the azimuth of the incident position is denoted by , and the interior angle between the ray and the direction of the z axis is denoted by . Then a specific incident ray at point A can be expressed in terms of ( , , ), also shown in Supplementary Figure 3b (iii). With traditional GRIN lens ray tracing 16 , the ray path (which represents the radius at distance ) of the specific incident beam can be represented by: where is the propagation length along the axis and is a constant determined by the manufacturing process that is directly related to the refractive index distribution, π √ is the period of the sinusoidal ray trace and the amplitude is ′ + . The refractive indices seen by the rays and rays are denoted by and respectively, which are functions of the radius . We then could express the refractive index of a GRIN lens in series form to approximate the real distributions: ( ) = (0) + + + + ⋯ + , ( = 1,2,3, ⋯ ) ( ) = (0) + + + + ⋯ + , ( = 1,2,3, ⋯ ) Where (0) and (0) are the refractive indices of the rays and rays at the centre (note that (0) = (0) ), and , , , ⋯ and , , ⋯ are constants determined by the manufacturing process. The effective refractive index ( , ) experienced by the ray at the local coordinates ( , , ) inside the GRIN lens is represented by: where is the interior angle between the wave normal and the extraordinary axis (blue double-sided arrow in Supplementary   Figure 3b (iv)), which is the complementary angle to . Along the sinusoidal ray trace, there will be an accumulated phase difference between the rays and rays when the beam reaches the back surface of the GRIN lens. So, there would be a different overall dependent on the traced ray C.
It should be noted here that there would also exist a minor beam split between rays and rays due to the birefringence properties of the GRIN lens, with which we are able to modulate the axial resolution of the GRIN lens imaging system as validated in Fig. 3b and 3c in the main article (the point spread function measurement process is adopted with reference to Ref [15]). This splitting is due to the differing sinusoidal optical path lengths of the rays (Ce and Co respectively). The parameter A in Eq. (9), which is related to the amplitude (and period) of the sinusoidal ray trace, is also different for o rays and e rays (in practice, the parameter should be represented by Ao and Ae, separately, which are determined by the corresponding refractive index profile, Eq. (10)). As the collection of o rays and e rays inside GRIN lens are associated with radial and azimuthal linear polarisation eigenmodes, the corresponding linear polarized light fields play a special role in all GRIN lens based systems. Based upon this understanding, we aim to calculate the optical path length difference (OPLD) and the phase difference determined by the retardance σ: Here ( ) and ( ) are the local refractive indices experienced by rays and rays, as a function of distance along the sinusoidal optical path, as an arc length integral from the original point on the front surface to the exit point on the back surface. Since the birefringence in the GRIN lens is very small, when calculating the accumulated retardance, we can make the approximation that the path of the rays and the rays is the same (as shown in Eq. (12)). We use parallel incident light ( = 0 at z = 0) throughout this work, and as the refractive index and birefringence profiles are rotationally symmetric, if we define the wavelength of the incident beam as , Eq. (13) would allow us to calculate the corresponding .
There are several points that should be considered when estimating the uncertainties of the beam split for different GRIN systems. In this validation work, based upon our polarization measurements above, we estimated that the GRIN lenses we used had maximum birefringence level around 10 -5 , which is in the range of plausible values for the lithium ion-exchange process used in their manufacture. We tested three different types of pitch 2 GRIN lenses, which, due to their manufacturing processes, had different birefringence properties and lengths. We found that they exhibited variations in beam splitting ranging between ~1.5μm and ~3.0μm. These variations are due to the different optical properties of each lens. We have assumed here that ( ) and ( ) follow a parabolic approximation, which one could likely refine for more accurate modelling of the phenomena.

Supplementary Note 3: Experimental set-up for polarization field measurements
Supplementary Figure 4 demonstrates the set-up for polarization field measurements. Here we used a LED light source (633nm, 3mW, Δ =20nm), a polarizer (P1) and a quarter waveplate (QWP1) to generate a light field of uniform polarization that was incident on the GRIN lens cascade. We used a quarter waveplate (QWP2) and a polarizer (P2) to measure the Stokes vectors of the light field by rotating QWP2 to four different angles [17][18][19] , following the process according to Ref [17][18][19]36].
The principal equations for calculation of the Stokes vector light field are shown in Eq. (14) and Eq. (15).
where is the Stokes vector of the incident light field, and are MMs of the corresponding polarizer and waveplate.
is the corresponding output Stokes vector for the fast axis orientation state of the quarter waveplate.
is the MM of the quarter waveplate for the fast axis orientation. is a n×4 matrix known as the instrument matrix 19 , which is derived from • . is the intensity information recorded by the camera.
Supplementary Figure 4. Set-up for the polarization field measurement. LED: source; P1, P2: polarizer; QWP1, QWP2: quarter waveplate; Camera: detector. The rest of the structure is the GRIN lens cascade, which combined different kinds of GRIN lenses with interstitial optical elements. Note that imaging optics have been omitted for clarity.

Supplementary Note 4: Theoretical and experimental validation for polarization light field characterization
Vector beams (VB) have attracted great interest for a range of applications that take advantage of their structured polarization [20][21][22] . They can be used to control the properties of a beam focus in microscopy 23 , for the demonstration of Möbius band-like topologies 24 , for electron acceleration 25 and material processing 26 Supplementary Figure 6 shows how different GRIN cascades can be used to generate a range of full Poincaré beams.
Supplementary Figure 6a illustrates a type of GRIN lens (coloured red) that has a larger magnitude retardance profile and is almost equivalent to a sequence of three of the first type GRIN lens (shown in yellow, as also used in main article Fig. 2).
Supplementary Figure 6b shows the generated beams using right hand circular polarization incident light.
Based on these results in Supplementary Figure 6b, we can make several observations: 1) The GRIN lens cascades are able to generate high radial order VBs due to the nonlinear variation of retardance along the radial direction ( Supplementary Figure 6b  structures (from (i) to (iv)) which use identical lenses (yellow one) or different lenses (red one). These cascades were also used in Fig. 2 in the main article.
These beams were generated using right hand circular polarized incident light.

Supplementary Note 5: Set-up for interferometry
We employed commonly used interferometric methods 35,36 to characterize OAM generation. Supplementary Figure 7 illustrates the Mach-Zehnder interferometer, in which we used a He-Ne laser (633nm, 2mW, with a Gaussian intensity distribution), a polarizer (P1) and a quarter waveplate (QWP1) to generate a uniform polarization incident light field at the input to the GRIN lens cascade. The light reflected by the first beam splitter (BS1) into the reference arm, passed off a silver mirror (M1) then was modulated by a half waveplate (HWP1). The beam was expanded then reflected again by another silver mirror (M2) before being combined with the beam from the experimental arm by the second beam splitter (BS2). Finally, the second quarter waveplate (QWP2) and the polarizer (P2) filtered the polarization state of the beam before measurement of the interferogram at the camera.

Supplementary Note 6: Theoretical and experimental basis of OAM component analysis
A light beam possessing a helicoidal phase-front, described by the structure of e , carries OAM of ћ (where is an integer and ћ is the reduced Planck's constant) [37][38][39][40] . Since this revelation in 1992 37 , various further applications of these beams have been developed, such as optical manipulation [41][42][43][44] , optical communications [45][46][47] , and in quantum and nano-optics [48][49][50][51] where ( ) is the amplitude distribution. When considering OAM, we are only concerned about the phase term (2 + ) in Eq. (22), where ( ) is an initial phase delay determined by the retardance profile ( ) (see Supplementary Note 2). The exponent i (2 ) shows that the analyzed beam exhibits two units of OAM, which reveals that the GRIN lens can be used as a spin-to-orbital angular momentum convertor. Corresponding results can be found in Fig. 2a (iii) to (v), Supplementary   Figures 8a and 8c.
We now examine the properties of the generated vector = • in another two specific eigenpolarization bases: From calculation, [ ] can be written as exponential form: = e + e + e ( ) + e ( ) , = e + e + e + e .
In Supplementary Figure 8 below, the corresponding phase profiles, intensity profiles, as well as interference patterns using analyzers described by [1, 0] (both in experiment and simulation) can be found in Supplementary Figures 8d (iii) to 8h (iii) and 8d (vii) to 8h (vii). We also modified the PSA to obtain a wider variety of wave fronts for illustration; the sequence is shown in Supplementary Figure 8b  Similar to the polarization fields shown in the main article (Fig. 2), some demonstrations of the complex phase modulation abilities of higher order GRIN lens cascades are also shown in Supplementary Figures 8i to 8l. It can be seen that higher order cascades are able to generate different complex phase patterns that also feature higher order OAM components by choosing the input SOP and the analysis polarization state. An example can been seen in Supplementary Figure 8j, We calculate 1 ′′ and 2 ′′ by through the same process used before, as shown in Eq. (27) and (28). In the following Supplementary Figures 9a and 9b we demonstrate phase profiles, intensity distributions and interference patterns (where the object beam is interfered with the same SOP light).
The simulations in this section were based on the above theoretical equations. It should be noted that other combinations of incident polarizations or analysis eigenvectors, or even higher order GRIN lens cascades cases (combined with Eq. (18)), can be easily calculated through the same process above.
We have demonstrated here that the GRIN lens cascades have the potential to generate complex structured phase profiles under different input/output eigenpolarizations. This novel structure for phase modulation (including the generation of OAM), derives from the special birefringence profile of GRIN optics. To our knowledge, this has not yet been reported. It might pave the way for new methods of complex beam engineering and benefit corresponding applications such as OAM related quantum and nano-optics.
Supplementary Figure 9. Phase profiles, interference patterns, intensity distributions generated by the single GRIN lens. There are many ways to specify the illumination polarizations and analyser configurations for MM polarimetry 1, 3-5, 7-11, 13, 18, 52, 53 . To find optimum measurement configurations, one can use the condition number (CN) of the system as an optimisation criterion for both the PSG and PSA. The CN can also be related to the volume of the inner tetrahedroid constructed by the SOPs inside the Poincaré sphere (shown in Fig. 4b (ii) of the main article) 18,52 .
Based upon this, we chose a PSG that consisted of four sectors of differently oriented quarter waveplates, as shown in Fig. 4b (i) in the main article, where the fast axes of waveplates 1 to 4 were oriented at 45°, 15.9°, 74.1°, and -45° according to the optimization shown in Ref [52]. Similar CN analysis has shown that the 132° phase retardance with four different fast axis directions is optimal for the PSA 18, 52 , providing a CN of 1.732. As shown in Fig. 4c Where the intensity information recorded by the camera for each modulation and detection channel are defined as a vector: and is defined by: which is a 4  4m matrix that consists of columns that are individual Stokes vectors from each combination of PSG sectors and camera pixels. Then the MM of the sample can be calculated as: where in −1 is the pseudo inverse matrix of in , and −1 is the pseudo inverse matrix of .

Supplementary Note 9: Validation of feasibility using standard samples
In order to test the capabilities of the MM polarimeter, we used as samples four polarizers with different orientations (0°, 90°, 45°, -45°), which are illustrated as P2 to P5 in the experimental setup in Supplementary Figure 10a. We took measurements of each polarizer in a single-shot and compared the derived MM elements with the ground truth MMs (Supplementary Figures   10b (i) to 10b (iv)). The maximum errors of the derived MM elements were smaller than 6.93%. The GRIN lens cascade based polarimeter has the following advantages: 1) it can perform single-shot MM measurement in a robust and stable manner; 2) single-shot measurements enable precise measurement of moving/changing objects; 3) it has the potential to be miniaturised into an integrated instrument, especially a fibre-based probe with scanning detection for clinical applications.  We measured 10 points per sample using our new MM polarimeter and a conventional MM microscope (overall 20 samples), as ground truth for quantitative comparison. In each sample, we chose the points by using MM polarimeter with the diameter of the field of view (FOV) at 0.19mm (Femto Technology Co. Ltd., G-B151157-S1483). The FOV of the MM microscope (calibration precision < 0.3%) was around four times larger than that of the polarimeter, so we took measurements from chosen sub-areas within same FOV of MM polarimeter to be processed. Then we calculated the mean value of the retardance across these areas to set as ground truth for further comparison through the same quantitative comparison process used to differentiate the breast cancerous stages by polarization parameters 53 . Example MMs, as well as corresponding MMPD parameters, from healthy or cancerous tissue are illustrated alongside quantitative statistic histograms in Supplementary   Figure 11.
For simplicity of measurement by the MM polarimeter we make the assumption that measurement of the average SOP should be sufficient for discrimination purposes. We first defocused the GRIN lens image, in order to implement optical averaging of the polarization effects of spatially varying structures. Then we began with a chosen region with optimal CN (related to the 132° retardance 18 ) to calculate the MM for demonstration of feasibility. If strong non-uniform area occurred such as near to the boundary between tissue types, we chose alternative rings (corresponding to other retardance values within the effective CN region on the PSA 52 ) or applied another single-shot to reconstruct the MM. In principle, alternative region shapes can be chosen for optimization of the processes.
All the obtained data measured by the MM polarimeter were used to draw an overall 3-D dot distribution in Fig. 4f in the main article. Then, combined with the corresponding MM microscope data, Supplementary Table 1 was formed to show the data measured by two approaches as well as the corresponding P-value 53 , which shows the significant difference between the two classes of samples. Considering the projection on Y (sample number) -Z (retardance) plane in the main article (Fig. 4f) and the P-values in Supplementary Table 1 of either individual samples or the overall combination, it can be seen that the polarimeter is able to distinguish efficiently the healthy and cancerous tissues. The difference between two kinds of samples is obvious across the tested samples; this also confirms the robustness of the polarimeter. Compared to the data measured by the MM microscope, larger standard deviations from the data measured by the polarimeter can be observed. The difference might result from the inhomogeneous sample structure in the effective FOV. To make the process more efficient and precise, further detailed error analysis and corresponding optimization will be the subject of further work. *the P-value is obtained by the significance test method 53 (P<=0.05 is significant), which shows the significance of the difference between the two sets of data.