Sorting OAM modes with metasurfaces based on raytracing improved optical coordinate transformation

Optical coordinate transformation (OCT) has attracted widespread attention in the field of orbital angular momentum (OAM) (de)multiplexing or manipulation, but the performance of OCT would suffer from its distortion. In this paper, we quantitatively analyze the distortion of OCT from the perspective of ray optics, and explain its rationality to work under non-normal incident light. For the special case of log-polar coordinate transformation (LPCT), we use a raytracing assisted optimization scheme to improve its distortion, which is related to a Zernike polynomial based phase compensation. After raytracing optimization, the root mean square error (RMSE) of the focused rays is reduced to 1/5 of the original value and the physical optic simulation also shows great improvement. In the experiment, we use three phase masks which are realized by metasurfaces, the measured results show well consistency with the simulation. Results in this paper have great potential to improve the performance of OCT related applications.


Introduction
In recent years, the ever-growing demands for information bandwidth has posed a huge challenge to the current optical communication system. Since the time, frequency, wavelength and polarization dimensions of light have been deeply researched before, exploiting the space dimension of light, also known as space division multiplexing (SDM), is considered to be the next-generation communication technology [1,2]. Among SDM, one of the most representative subsets is orbital angular momentum (OAM) multiplexing, which have given rise to especial research interest in mode division multiplexing (MDM). Generally, light beam with spiral phase wavefront factor , where is the azimuthal coordinate, is the topological charge (TC), carry ℏ OAM per photon [3]. In principle, is an unbounded integer, OAM beams with different are orthogonal and thus have great application prospects in optical communication [4]. Apart from communication, OAM beams have also been applied to fields such as optical manipulation [5], optical tweezer [6], imaging [7], and quantum information [8]. In all of these applications, is a determined parameter since it is directly related to the quantity of OAM. Owing to its great importance, measuring or demultiplexing different OAM beams have attracted extensive researches [9][10][11][12][13][14][15].
Although there have been many studies on OAM measurement, simultaneously demultiplexing multiple OAM beams without shrinkage of efficiency is still challenging. To cope with this problem, in 2010, Berkhout et al. proposed a scheme which is based on optical coordinate transformation (OCT). In their way, annular OAM beam with spiral phase is firstly transformed into rectangular beam with plane wave phase and then can be sorted by a focused lens, this method also known as log-polar coordinate transformation (LPCT) [16]. Due to the simple implementation and unified efficiency, the LPCT method has received widespread concerns [17][18][19][20][21][22]. Inspired by the LPCT, some other versions of OCT for OAM manipulation have also been proposed, such as circular-sector transform [23] and spiral transform [24]. Though the OCT has its great advantage in sorting OAM, its performance could suffer from severe distortion when used for higher-order OAM state, this is because the beam need to assumed to be normally incident while the higher-order OAM beam has large skew angle [16]. To relieve this problem, Lavery et al. have analyzed the distortion and used two large phase masks to achieve up to 28 orders of OAM measurement in 632nm [19]. However, though we know the distortion could be introduced by non-normal incident light rays, there are still no quantitively analysis of the distortion and no corresponding solutions to improve the performance of the OAM sorter that work under specific configuration.
In this work, through analyzing the ray's trajectory among the coordinate transformation process, we will explain the rationality to use the OCT beyond normal incident assumption. It should be noted that the ray deviated from normal incidence will go through incorrect deflection in the second phase plate. Then we can further provide a quantitative result of the distortion in OCT, and for the case of LPCT, we give an upper bound of that can be measured under certain distortion. In addition, to alleviate the performance degradation caused by the distortion, we present a scheme to improve the distortion by raytracing assisted optimization, which is related to Zernike polynomial based phase compensation. After the optimization, we use three phase masks realized by metasurfaces to build the experimental setup and implement the measurement.

Theoretical analysis of distortion
Generally, the OCT maps the points of two separate planes, let ( , ) denote the first plane and ( , ) denote the second. Assuming light is normally incident, to achieve the map, from the view of ray optics and the paraxial approximation condition, a phase mask should be put on the first plane and its gradient satisfied 22 , where is the wavelength, and is the distance between the two planes. When the ray arrives at the second plane, similarly, we need a phase mask to compensate the ray's deflection, and its gradient satisfied In general, given a reasonable relationship between ( , ) and ( , ), we can figure out the corresponding and from above equations to achieve the transformation. As shown in Fig.  1, then a ray normally incident on the ( , ) in the first mask, will pass through the device and arrives at the ( , ), and normally incident out from the second mask. However, if the ray is non-normally incident, its trace will deviate from the position ( , ), and also experience unexpected phase gradients. In this situation, further discussion is needed. For simplicity, assuming the incident ray can be described by phase distribution , then the ray's deviation in the second mask caused by its phase gradient can be described as This means that if the incident light is non-normal incidence or ∇ ≠ 0, it will cause the distortion of position of OCT. Since OAM beam has skew angle, to deal with this problem for better performance, | ⃗ | ≪ is supposed in reported works [17,24] and thus can cause a restriction on the of OAM beam with given size . In above, the position after transformation is considered within the scope of ray optics. However, for the application of OAM (de)multiplexing or manipulation, the distortion of position after OCT is not the only problem, because the phase transformation is the main purpose. To better evaluate the distortion, ray's direction after OCT is need to be considered. Due to the position distortion, the ray will obtain a corresponding phase gradient change in the second plane, and by using Eq. (2) the gradient change can be described as where ∆ = ( + , + ) − ( , ), ∆ = ( + , + ) − ( , ), since the ray's deflection caused by the first mask can be offset by the second mask, its ultimately phase gradient accumulated by the transformation process can be expressed as Eq. (5) shows the phase gradient that represents the ray's direction, to explore the ray's property, we perform Taylor expansion on ⃗ and the gradient can be written as where ∇′ is arithmetic of ( , ) coordinate, ∆ ⃗ is higher order residuals. From Eq. (2) and using ⁄ = ⁄ , the ray's phase gradient can be further written as '. G P G =  +  (7) From the right hand of Eq. (7), we can see that the input ray with phase gradient ∇ has transformed into the output ray with gradient ∇′ after OCT, it is also mean that the phase distribution in the first plane has transformed into the second plane. In this way, the rationality of OCT to transform the phase distribution is justified by the first-order Taylor approximation. Meanwhile, the residual term in the equation is actually a distortion that cause the diffusion of the ray and can be calculated from Eq. (5) and Eq. (7). For convenience, this distortion can also be approximated by the second order Taylor expansion In most situations, the coordinate transform is conformal, or more specifically, it is a conjugate analytic function ( ) where is an analytic function, then using the distortion can also be expressed as From Eq. (10), we can conveniently calculate the distortion of various OCT. From Eq. (3) we can see the position deviation is proportional to | ⃗ |, while here the phase gradient distortion |∆ ⃗ | is proportional to | ⃗ | 2 . Next, we start to consider the special case of LPCT that have widely used for OAM sorting. In this case, = − ln( / ), = , where , are positive parameters of the transformation, then we can derive the corresponding relationship (11) In addition, for an OAM beam with given radius and phase wavefront = , from Eq. (3) the deviation caused by its skew angle can be described as Then from Eq. (10), the distortion can be expressed as Eq. (13) is the distortion of LPCT, if we postulate the distortion is smaller than the difference between adjacent OAM modes, because otherwise the crosstalk will become too large, that's mean |∆ ⃗ | (1/ ) ⁄ ≤ 1, and we can get This is the upper bound of the OAM states can be sorted by the LPCT under given distortion. For other type coordinate transformations such as spiral transformation [24], through same treatment, similar constraint can also be obtained.

Raytracing assisted distortion improvement
In order to improve the distortion of LPCT derived in above section, we introduce the raytracing to simulate the process. Since the raytracing is a common tool for geometrical optics design, it should also be suitable for analyzing OCT. Here, to implement the LPCT, we use three phase masks to work as transformer, compensator, and a lens respectively to separate or sort the input OAM beams. These three phase masks are shown in Fig. 2(a-c). In addition, the input OAM beam used here is LG beam, of which beam size is varied according to and follow = √| |, where is the waist radius of the gauss mode. In this situation, to limit the interference on adjacent order as derived in above section, the range of TC should satisfy The raytracing then is done according to Eq. (1-2) where paraxial condition is used. Fig.  2(d) shows the raytracing schematic diagram of sorting multiple OAM beams, different OAM beam has corresponding skew angle determined by its order and radius, and the rays will go through the three masks and convergence at the focal plane. As is shown in Fig. 2(e), the rays of OAM beams with different order , where −10 ≤ ≤10, have different radius. The waist of the LG beam is set to be 0. 18 , and the parameters of the LPCT respectively set as, = 1.55 , = 1.5/2 , = 0.45 , and = 15 . For simplicity, the focal length of the convergence lens is also set to be . Then from Eq. (15) we can find the critical TC is about = 8. Accordingly, Fig. 2(f) shows the results of focused rays based on the installation, it can be seen that the ray will go to adjacent focus when the order is larger than 8. This ray optics simulation result confirmed our theoretical analysis and shows good accuracy to approximate the distortion of LPCT using second order Taylor expansion. Next, to overcome the huge distortion of LPCT resulted by the skew angle of input OAM beam, we further use the raytracing to optimize the phase masks. To implement the optimization, we define a Zernike polynomial based phase compensation. For every mask, the first 45 order Zernike polynomials are added to express the phase compensation, fewer terms of polynomials may lead to worse results while more terms would increase the difficulty of optimization. The modified phase mask thus can be described as follow where is the original phase, is the -order Zernike polynomial, is the corresponding coefficient, and ' is the modified phase. When a ray goes through the mask, its traces will be regulated by extra phase denoted by Zernike polynomial, and the coefficients thus can be used as optimizing variables. Besides, we define the root mean square error (RMSE) as the extent of distortion for the optimization, this can be denoted by  (17) where needs to be minimized, is the number of rays, ( , ) and ( ,̃) are the ray's arrival position and ideal destination in the focal plane. After set up this raytracing model, we use an evolutionary algorithm named covariance matrix adaptive moment estimation (CMA-ES) to regulate the optimizing parameters or the Zernike coefficients [25], which is a commonly used global optimization algorithm. Other methods such as genetic algorithm (GA) [26], or particle swarm optimization (PSO) [27] may be also valid here. Figs. 3(a-c) show the phase masks after Zernike polynomial based phase compensation. Fig.  3(d) shows the corresponding convergence of rays in focal plane. There is a great deal of diminution of distortion especially for the higher-order OAM, and the RMSE have dropped from original value of 15.87 to 3.27 while the spacing of adjacent TC is 15.50 . In Figs. 3(e-f) and Figs. 3(g-h), we show the physical optics based simulation results before and after optimization. Fig. 3(e) and 3(f) show the intensity distribution of multiple OAM beams in the focal plane and corresponding normalized intensity distribution along the central axis. From these focused spots we can find that with the increasing , there are also increasing distortion. The critical order is = 8 when the distortion goes to heavily influence the adjacent spots, which is in good agreement with the theoretical prediction. Correspondingly, the simulation results after optimization are shown in Fig. 3(g) and 3(h). From the figures we can see a huge improvement of the sorting performance. The spots in Fig .3(g) have roughly uniform size and less overlap with adjacent spots. Meanwhile, the intensity distribution of the sorted OAM beams along the central axis in Fig. 3(h) also indicate less interference to other beams. These simulation results signify that the distortion improvement is effective.

Device preparation and characterization
In order to prepare above optimized phase masks for experiment, we use three metasurfaces to provide the corresponding phase modulation. Here the function of metasurfaces is mainly concentrate on 2π phase coverage and high transmittance, a most favorable selection is the dielectric metasurface that have demonstrated in reported researches [28][29][30]. As shown in Fig.  4(a), the metaunit of our dielectric metasurface is composed of amorphous silicon (a-Si) cylinder and SiO2 substrate, the metaunit thus is polarization insensitive. In addition, to prevent the cylinder from collapsing in the fabrication, 20nm a-Si is remained in the substrate. Under the given structure, to maximize the transmittance while cover 2π phase, we have done a parameter sweep to determine the metaunit and respectively set = 760 and ℎ = 890 according to the sweeping results. The phase and transmittance response of the metaunit are shown in Fig. 4(b). We select 96 structures with phase difference 2 /96 through interpolation and have marked them with asterisk in the Fig. 4(b). Based on the metaunit, coverage of 2 phase can be achieved and the average transmittance of these structures is 90.77%. After the metaunit design, we generate the layout to meet the required phase masks in above section and fabricate the device accordingly. The fabricated process can be briefly described as follow. Firstly, a layer of 890-nm-thick a-Si is evaporated in the SiO2 substrate using chemical vapor deposition (CVD). After pretreatment of the sample, the AR-P6200.13 (ARP) resist is spin coated in the a-Si film. Before electron beam lithography (EBL), a layer of Al is deposit in the ARP resist. The device then is put into the EBL device to write the mask. After exposure, we remove the Al layer and develop the resist to get the patterned mask. Afterwards, we use inductive coupled plasma (ICP) to etch the a-Si layer. Finally, the ARP resist is removed and the device is fabricated. In Figs. 4(c-e), we show the fabricated metasurfaces that are captured in microscope, which are consistent with Figs. 2(a-c). Figs. 4(f-g) show the partial fabrication detail in scanning electron microscope (SEM). From these images we can see the device is well fabricated and the pillars have been etched appropriately. Based on the metasurfaces described above, we built the experiment measurement setup as shown in Fig. 5. A laser with 1550nm wavelength is used as the light source and the generated beam is collimated into SLM, where we use a PC to regulate its polarization. After the reflection and modulation of the SLM, the OAM-carrying beam then is contracted by a pair of lenses L1 and L2 to match the size of the metasurfaces. To diminish the unexpected diffraction order of SLM, an iris is placed between the two lenses. The OAM beam then arrive at the metasurface MS1 and undergo the modulation of MS1-3. Finally, since the focal spot is small, we use lenses L3 and L4 to expand it and the CCD is used to detect the results. At last, we carry out the experiment, the sorting results of OAM beams with −16 ≤ ≤ 16 are shown in Fig. 6. As comparison, Figs. 6(a-b) show the simulation results before and after optimization. It can be seen that the simulation results have severe distortion before optimization when l≥ 8, but after optimization, the divergence of focal spot are improved significantly. As shown in Fig. 6(c), the experimental detection of OAM beams up to ±16 order is achieved, and the distortion of higher-order OAM beams also have been significantly suppressed. Further, in order to quantify the effect of distortion improvement, the detected transfer matrix of simulation and experiment are shown in Fig. 7. Since the focal spot is slightly inclined, to better evaluate the improvement of the spot divergence, the relative power of the detected transfer matrix is determined by mean intensity along the vertical axis in focal plane. In this way, from Fig. 7(a) we can see that the crosstalk is very high before optimization, and the higher -order OAM beams become gradually indistinguishable with adjacent beams when l≥ 8. This situation can be improved by raytracing assisted optimization, as shown in Figs. 7(b-c), both the simulation and experimental results show well distinction even for l= ±16. Though there are some differences between the simulation and experiment results, which may be due to the manufacturing errors and the imperfect experimental setup. On the whole, the performance of LPCT has been improved greatly, and the experimental results can well support the reliability of simulation and indicate a good effectiveness of the raytracing assisted distortion improvement.

Conclusions
In conclusion, we have provided a quantitative analysis of the distortion in OCT, and also clarified the rationality of optical coordinate transformation for non-normally incident beams in the view of geometrical optics. For the case of LPCT, we give an upper bound of l that can be measured under certain distortion. As the OCT can only work well in the first order Taylor expansion approximation, the sorting distortion will increases as the TC of OAM beam increases. Therefore, to alleviate this problem, we further provide a method based on the raytracing and Zernike polynomial based phase compensation. After raytracing optimization, the RMSE of the focused rays can be reduced to 1/5 of the original and the physical optic simulation also shows obviously distortion improvement. At last, to confirm the simulation results, we implement the experiment using three fabricated metasurfaces. The experimental results are in good agreement with the simulation and show obvious performance improvement after raytracing optimization. Results in this paper may be also help to further improve other OCT related applications and have potential to extend by adding more masks.