Numerical Simulation of Integrated Generation and Shaping of Airy and Bessel Vortex Beams Based on All-Dielectric Metasurface

Integrating multiple independent functions into a single optical component is one of the most important topics in research on photoelectric systems. In this paper, we propose a multifunctional all-dielectric metasurface that can achieve a variety of non-diffractive beams depending on the polarization state of the incident light. Using the anisotropic TiO2 rectangular column as the unit structure, the three functions of generating polygonal Bessel vortex beams under left-handed circularly polarized incidence, Airy vortex beams under right-handed circularly polarized incidence and polygonal Airy vortex-like beams under linearly polarized incidence are realized. In addition, the number of polygonal beam sides and the position of focal plane can be adjusted. The device could facilitate further developments in scaling complex integrated optical systems and fabricating efficient multifunctional components.


Introduction
In 1992, Allen et al. [1] discovered for the first time that photons can carry orbital angular momentum (OAM) in addition to general spin angular momentum (SAM) [2], which laid the foundation for the study of vortex beams. Optical vortex beam (OVB) is a special beam carrying an OAM of lh per photon. Due to its singularity effect, OVB has a helical phase exp(ilϕ) in the wavefront structure, where l is the topological charge (TC) number which can take any integer value, ϕ denotes the azimuthal angle andh is the Planck constant. Since then, OVB has been widely used in quantum information processing [3,4], quantum entanglement [5], optical micromanipulation [6,7], optical communication [8][9][10] and nonlinear optics [11,12]. Hence, with the increasing demand for practical applications, OVB needs to be further studied: from generating [13][14][15] and measuring [16] to shaping OVB [17]. Especially in recent years, shaping OVB has a research hot spot. Current common shaping methods such as multi-OVB interaction [18][19][20] and optical oscillatory elements [21] require complex holography or high-precision diffractive elements. Until 2019, a new type of phase structure called higher-order cross-phase (HOCP) was introduced into the modulation of the beam. The function of HOCP is to introduce astigmatism into OVB in phase, so that the optical field mode changes during the propagation process, thereby changing the shape of OVB [22]. Compared with previous methods, it has higher integration and simpler implementation, which opens up a new horizon for shaping OVB.
However, in practical applications, a single OVB will be affected by many clutters in the propagation, which will deteriorate the performance. In order to further stabilize and improve the performance of OVB, a combination of OVB with other beams is a solution. In 1979, Berry and Balazs [23] first theoretically proved that the Schrödinger equation has a solution in the form of an Airy function and predicted a non-diffractive Airy wave packet. The optical

Design of Multifunctional Metasurface
In this section, we will describe the physical mechanism that achieves the polarization sensitivity difference in the x and y direction in the transmission mode. We consider an arbitrary anisotropic nanopillar unit placed parallel to the z-axis, with the long and short axes along the x and y directions, respectively. When the unit is excited by the incident light with different polarization states, it produces responses with different phases and amplitudes along the x and y directions. Under LP light incidence, the Jones matrix of the transmitted field for an arbitrarily anisotropic unit rotated by an angle around the z-axis in the xoy plane is expressed as follows: where Tx, Ty, φx and φy represent the transmission amplitudes and phases along the x and y axes. θ is the rotation angle of the unit relative to the reference coordinate system. R represents a 2 × 2 rotation matrix. Here, we assume that the transmission amplitude of the unit remains unity amplitude (Tx = Ty = 1) with a π phase difference (φx = φy + π) for the x and y polarization cases.
The independent phase modulation of LCP and RCP on the metasurface can realize the generation of φ1(x, y) and φ2(x, y) two different phase distribution beams. This means that when the RCP is incident, the metasurface will only be modulated by φ1(x, y). Similarly, the metasurface will only generate a beam modulated by φ2(x, y) when the LCP is incident. This modulation mechanism on metasurface can be represented by a matrix J(x, y): After the matrix transformation in Equations (2) and (3), the matrix J(x, y) can be expressed as:

Design of Multifunctional Metasurface
In this section, we will describe the physical mechanism that achieves the polarization sensitivity difference in the x and y direction in the transmission mode. We consider an arbitrary anisotropic nanopillar unit placed parallel to the z-axis, with the long and short axes along the x and y directions, respectively. When the unit is excited by the incident light with different polarization states, it produces responses with different phases and amplitudes along the x and y directions. Under LP light incidence, the Jones matrix of the transmitted field for an arbitrarily anisotropic unit rotated by an angle around the z-axis in the xoy plane is expressed as follows: where T x , T y , ϕ x and ϕ y represent the transmission amplitudes and phases along the x and y axes. θ is the rotation angle of the unit relative to the reference coordinate system. R represents a 2 × 2 rotation matrix. Here, we assume that the transmission amplitude of the unit remains unity amplitude (T x = T y = 1) with a π phase difference (ϕ x = ϕ y + π) for the x and y polarization cases.
The independent phase modulation of LCP and RCP on the metasurface can realize the generation of ϕ 1 (x, y) and ϕ 2 (x, y) two different phase distribution beams. This means that when the RCP is incident, the metasurface will only be modulated by ϕ 1 (x, y). Similarly, the metasurface will only generate a beam modulated by ϕ 2 (x, y) when the LCP is incident. This modulation mechanism on metasurface can be represented by a matrix J(x, y): The matrix form of |LCP corresponds to 1 i , and |RCP corresponds to 1 −i . After the matrix transformation in Equations (2) and (3), the matrix J(x, y) can be expressed as: J(x, y) = e iϕ 1 (x,y) e iϕ 2 (x,y) −ie iϕ 1 (x,y) ie iϕ 2 (x,y) Then, by making the transmission matrix T(x, y) of Equation (1) equivalent to the modulation matrix J(x, y) of Equation (4), the phase distribution along the x and y directions and the unit rotation angle can be calculated, respectively.
The arrangement of units in the metasurface is determined by Equations (5)- (7). The phase design is mainly determined by the combination of the Pancharatnam-Berry (PB) phase adjustment based on Equation (5) and the propagation phase adjustment based on Equation (7). The PB phase is based on the adjustment of the phase by rotating the direction of the half-wave plate under circularly polarized light, which can be explained by the Poincaré sphere [51]. As shown in Figure 2, the north and south levels of the Poincaré sphere, respectively, represent LCP and RCP, and each point on the equator line corresponds to a different linear polarization state, while the entire spherical surface The other points represent elliptical polarization states for different polarization states. When a point on the sphere rotates an angle Ψ with the center of the sphere as the center, the phase change shown on the sphere is twice the rotation angle, which is the threedimensional interpretation of the geometric phase. Therefore, when the rotation angle of nanopillars varies from 0 to π, the PB phase of 0 to 2π can be obtained. Since this phase is only related to the orientation and size of units, and since the different polarization states of the incident light will produce different phase responses, this feature can be used to perform differential encoding of left-handed and right-handed incidence, thereby realizing multifunctional metasurfaces.
Then, by making the transmission matrix T(x, y) of Equation (1) equivalent to the modulation matrix J(x, y) of Equation (4), the phase distribution along the x and y directions and the unit rotation angle can be calculated, respectively.
The arrangement of units in the metasurface is determined by Equations (5)- (7). The phase design is mainly determined by the combination of the Pancharatnam-Berry (PB) phase adjustment based on Equation (5) and the propagation phase adjustment based on Equation (7). The PB phase is based on the adjustment of the phase by rotating the direction of the half-wave plate under circularly polarized light, which can be explained by the Poincaré sphere [51]. As shown in Figure 2, the north and south levels of the Poincaré sphere, respectively, represent LCP and RCP, and each point on the equator line corresponds to a different linear polarization state, while the entire spherical surface The other points represent elliptical polarization states for different polarization states. When a point on the sphere rotates an angle Ψ with the center of the sphere as the center, the phase change shown on the sphere is twice the rotation angle, which is the three-dimensional interpretation of the geometric phase. Therefore, when the rotation angle of nanopillars varies from 0 to π, the PB phase of 0 to 2π can be obtained. Since this phase is only related to the orientation and size of units, and since the different polarization states of the incident light will produce different phase responses, this feature can be used to perform differential encoding of left-handed and right-handed incidence, thereby realizing multifunctional metasurfaces.

Design of the Unit Structure
In order to better realize the combination of the transmission phase and the PB phase, the design of the units in the metasurface needs to be adjusted by the rotation angle and size. The design of a unit is shown in Figure 3. Due to the high refractive index and relatively low loss of TiO2 dielectric material at visible wavelength, TiO2 semiconductor is

Design of the Unit Structure
In order to better realize the combination of the transmission phase and the PB phase, the design of the units in the metasurface needs to be adjusted by the rotation angle and size. The design of a unit is shown in Figure 3. Due to the high refractive index and relatively low loss of TiO 2 dielectric material at visible wavelength, TiO 2 semiconductor is selected as the rectangular unit material which placed on the SiO 2 substrate. The refractive index parameters (n, k) of both TiO 2 and SiO 2 materials are taken from [52]. To further discover the realization mechanism of phase, we performed numerical simulations using a three dimensional finite difference time domain (FDTD) solutions. We use scripts to import the above theoretical content into FDTD, and obtain corresponding results under 3D simulation to support the above theoretical derivation. The size of a mesh cell we set is 0.02 µm, and the simulation wavelength is set as λ = 532 nm. In the case where the simulation sets this value to the wavelength condition, the real part (n) of refractive indices of TiO 2 and SiO 2 are n TiO 2 = 2.17 and n SiO 2 = 1.48, and the imaginary part (k) of the refractive indexes are both approximately zero. Considering both integration and transmission effects in this design, we set H = 900 nm and P = 400 nm, and L and W were varied from 50 nm to 400 nm, respectively. Figure 3b shows the normalized transmittance and phase magnitude of the 11 units selected after the simulation. It can be clearly seen that the selected units satisfy the full coverage of the 2π range well, and has a constant phase response difference of approximately π for x and y polarization. We apply periodic boundary conditions in the x and y directions and implement a perfectly matched layer (PML) in the z direction. The unit can be thought as a rectangular waveguide with birefringent properties due to the asymmetry of its geometry. Therefore, the phase delay along the x and y axes can be achieved by changing L and W.
selected as the rectangular unit material which placed on the SiO2 substrate. The refractive index parameters (n, k) of both TiO2 and SiO2 materials are taken from [52]. To further discover the realization mechanism of phase, we performed numerical simulations using a three dimensional finite difference time domain (FDTD) solutions. We use scripts to import the above theoretical content into FDTD, and obtain corresponding results under 3D simulation to support the above theoretical derivation. The size of a mesh cell we set is 0.02 µm, and the simulation wavelength is set as λ = 532 nm. In the case where the simulation sets this value to the wavelength condition, the real part (n) of refractive indices of TiO2 and SiO2 are , and the imaginary part (k) of the refractive indexes are both approximately zero. Considering both integration and transmission effects in this design, we set H = 900 nm and P = 400 nm, and L and W were varied from 50 nm to 400 nm, respectively. Figure 3b shows the normalized transmittance and phase magnitude of the 11 units selected after the simulation. It can be clearly seen that the selected units satisfy the full coverage of the 2π range well, and has a constant phase response difference of approximately π for x and y polarization. We apply periodic boundary conditions in the x and y directions and implement a perfectly matched layer (PML) in the z direction. The unit can be thought as a rectangular waveguide with birefringent properties due to the asymmetry of its geometry. Therefore, the phase delay along the x and y axes can be achieved by changing L and W.  Figure 4a,c depict the relationship between the size (L, W) of units and the phase shift and transmission coefficient under x-LP light, respectively. It can also be seen that the phase retardation can well span the entire 2π range, and the transmittance is always at a high value of more than 80%. Therefore, the desired phase combination can be obtained by rationally choosing the unit size. In addition, the selected units are well satisfied that the phase difference are approximately π under x-LP and y-LP incidence, which can be clearly drawn from Figure 4b. The calculated phase shifts and transmission coefficients of individual TiO2 unit under x-LP and y-LP light incidence are presented in Figure 4d. In a word, the independent regulation of RCP and LCP can be well achieved by all units.  Figure 4a,c depict the relationship between the size (L, W) of units and the phase shift and transmission coefficient under x-LP light, respectively. It can also be seen that the phase retardation can well span the entire 2π range, and the transmittance is always at a high value of more than 80%. Therefore, the desired phase combination can be obtained by rationally choosing the unit size. In addition, the selected units are well satisfied that the phase difference are approximately π under x-LP and y-LP incidence, which can be clearly drawn from Figure 4b. The calculated phase shifts and transmission coefficients of individual TiO 2 unit under x-LP and y-LP light incidence are presented in Figure 4d. In a word, the independent regulation of RCP and LCP can be well achieved by all units.

Polygonal Bessel Vortex Beams
The phase φ1(x, y) of the cross phase (CP) in Cartesian coordinates (x, y) has the form [22]: where the parameter u controls the conversion rate, the azimuth factor θ characterizes the rotation angle of the transformed beam on the Fourier plane, p and q are arbitrary positive integers, and the value of p plus q is equal to the order of CP. When CP is low-order cross phase (LOCP), that is, both p and q have a value of 1, the number of topological charges can be determined. Further, when the sum of p and q exceeds 2, CP as the high-order cross phase (HOCP) can realize adjusting the beam to an arbitrary polygon. This paper mainly exploits the ability of the CP to shape the beams, so the HOCP mode is adopted. For computational convenience, we consider the case where the CP does not rotate (θ = 0). The phase φx(x, y) of adding HOCP to the Bessel vortex beam is: q p x y x u x y l y x y x (9) where λ represents the incident wavelength, β represents the angle between the transmitted wave and z-axis (set to 12.7°), l represents the topological charge, the parameter u controls the beam conversion rate and η is an arbitrary real number. It can be seen from Equation (9) that the magnitude of the value of u is affected by the order of CP. Since the designed structure is nanoscale, the value of u will be relatively large. As an example, in the following Figure 5, the generation process of the phase mask composed of the HOCP and the 2nd-order Bessel vortex beam is shown. Figure 5a1-a4 corresponds to the phase masks from the 3rd-order CP to the 6th-order CP, respectively, and Figure 5b1

Polygonal Bessel Vortex Beams
The phase ϕ 1 (x, y) of the cross phase (CP) in Cartesian coordinates (x, y) has the form [22]: where the parameter u controls the conversion rate, the azimuth factor θ characterizes the rotation angle of the transformed beam on the Fourier plane, p and q are arbitrary positive integers, and the value of p plus q is equal to the order of CP. When CP is low-order cross phase (LOCP), that is, both p and q have a value of 1, the number of topological charges can be determined. Further, when the sum of p and q exceeds 2, CP as the high-order cross phase (HOCP) can realize adjusting the beam to an arbitrary polygon. This paper mainly exploits the ability of the CP to shape the beams, so the HOCP mode is adopted. For computational convenience, we consider the case where the CP does not rotate (θ = 0). The phase ϕ x (x, y) of adding HOCP to the Bessel vortex beam is: where λ represents the incident wavelength, β represents the angle between the transmitted wave and z-axis (set to 12.7 • ), l represents the topological charge, the parameter u controls the beam conversion rate and η is an arbitrary real number. It can be seen from Equation (9) that the magnitude of the value of u is affected by the order of CP. Since the designed structure is nanoscale, the value of u will be relatively large. As an example, in the following Figure 5, the generation process of the phase mask composed of the HOCP and the 2ndorder Bessel vortex beam is shown. Figure 5a 1 -a 4 corresponds to the phase masks from the 3rd-order CP to the 6th-order CP, respectively, and Figure 5b 1 -b 4 shows the phase of the 2nd-order Bessel vortex beam after superimposing HOCP in Figure 5a 1 -a 4 . It can be clearly seen that the Bessel vortex beam can be shaped into any polygon after adding HOCP, where the number of polygon sides is equal to the order of HOCP. Another way of saying this is that when the order of HOCP is three, four, five and six, the resulting OVB is shaped into triangular, quadrilateral, pentagonal and hexagonal polygons. Based on the above theory, considering the influence of p and q values on beam shaping, the xoy plane electric field simulated at the focal position under RCP incidence is shown in Figure 6. and 4th-order Bessel vortex beams without adding HOCP, respectively. It can be observed that Bessel vortex beams of any order can be well generated. When HOCP is added, the electric field presents a polygonal ring of corresponding order, and the phase presents a polygonal distribution. Among them, some basic properties are the same as the traditional Bessel vortex light, for example: the radius of the beam will increase with the increase in the order of the Bessel vortex beam, and there will be secondary halos of various levels outside the primary halo. From Figures 6b 1 -b 5 and 6d 1 -d 5 , it can be clearly known that the adjustment of HOCP is based on the phase. By polygonal shaping the vortex part of the phase mask, the shape of the spatial beam can be changed, which can confirm the numerical results derived in Figure 5. Since the manipulated particle will move in the direction of the energy flow, as long as the values of p and q are controlled according to the desired result, the motion trajectory of the manipulated instance can be accurately determined. adding HOCP, where the number of polygon sides is equal to the order of HOCP. Another way of saying this is that when the order of HOCP is three, four, five and six, the resulting OVB is shaped into triangular, quadrilateral, pentagonal and hexagonal polygons. Based on the above theory, considering the influence of p and q values on beam shaping, the xoy plane electric field simulated at the focal position under RCP incidence is shown in Figure 6. Figure 6a1-d1 represents the electric field and phase distribution of the 2nd-order and 4th-order Bessel vortex beams without adding HOCP, respectively. It can be observed that Bessel vortex beams of any order can be well generated. When HOCP is added, the electric field presents a polygonal ring of corresponding order, and the phase presents a polygonal distribution. Among them, some basic properties are the same as the traditional Bessel vortex light, for example: the radius of the beam will increase with the increase in the order of the Bessel vortex beam, and there will be secondary halos of various levels outside the primary halo. From Figure 6b1-b5 and Figure 6d1-d5, it can be clearly known that the adjustment of HOCP is based on the phase. By polygonal shaping the vortex part of the phase mask, the shape of the spatial beam can be changed, which can confirm the numerical results derived in Figure 5. Since the manipulated particle will move in the direction of the energy flow, as long as the values of p and q are controlled according to the desired result, the motion trajectory of the manipulated instance can be accurately determined.

Airy Beams
Secondly, the phase of the Airy beam is: where τ = k/f, k = 2π/λ, λ represents the wavelength of the incident light, f is the focused wavelength, and m and n are arbitrary constants. To verify that the phase structure can realized the Airy beams, we assume the complex amplitude of the plane wave moving

Airy Beams
Secondly, the phase of the Airy beam is: where τ = k/f, k = 2π/λ, λ represents the wavelength of the incident light, f is the focused wavelength, and m and n are arbitrary constants. To verify that the phase structure can realized the Airy beams, we assume the complex amplitude of the plane wave moving along the z-axis: After Fraunhofer diffraction of light field, the expression for the outgoing electric field is: where z 1 is the distance from the observation screen to the diffraction screen. Setting z 1 = f and substituting Equations (10)- (12), the change in the electric field can be deduced: As can be seen from Equation (14), the light field of an Airy beam is affected by m and n. We set m = n = a, where a is an any constant. Corresponding to different a values under LCP incidence, the xoz and yoz electric field distributions and the intensity distributions at the focus position (the dotted line position in Figure 7a 1 -a 5 are shown in Figure 7. Figure 7c shows the intensity distribution of the light field along the x-axis at the focus. With the increase in the value of a, the focal distance will gradually approach the metasurface and the central intensity of the focal position will progressively increase. In other words, the Airy beam becomes more concentrated as the value of a increases. where z1 is the distance from the observation screen to the diffraction screen. Setting z1 = f and substituting Equations (10)- (12), the change in the electric field can be deduced: As can be seen from Equation (14), the light field of an Airy beam is affected by m and n. We set m = n = a, where a is an any constant. Corresponding to different a values under LCP incidence, the xoz and yoz electric field distributions and the intensity distributions at the focus position (the dotted line position in Figure 7a1-a5 are shown in Figure  7. Figure 7c shows the intensity distribution of the light field along the x-axis at the focus. With the increase in the value of a, the focal distance will gradually approach the metasurface and the central intensity of the focal position will progressively increase. In other words, the Airy beam becomes more concentrated as the value of a increases.

Airy Vortex-Like Beams
Next, according to the focus position of the Bessel vortex beam, we correspondingly select the Airy beam with a = 16 × 10 −7 . Firstly, considering the case of the general Bessel vortex beam, when LP light is incident, the light fields generated by the incident LCP and RCP will act simultaneously. We slightly adjusted the output intensity of the two beams in order to achieve better fusion under the incidence of LP light. Additionally, we adjust

Airy Vortex-like Beams
Next, according to the focus position of the Bessel vortex beam, we correspondingly select the Airy beam with a = 16 × 10 −7 . Firstly, considering the case of the general Bessel vortex beam, when LP light is incident, the light fields generated by the incident LCP and RCP will act simultaneously. We slightly adjusted the output intensity of the two beams in order to achieve better fusion under the incidence of LP light. Additionally, we adjust the two light fields to focus on the same plane, and there will be an interaction between the two light fields, resulting in an Airy vortex-like beam, which can be verified by Figure 8. The light field and phase distribution of the xoy cross-section at the focal position are shown in Figure 8b,c. It can be found that after the superposition of the light fields, the central spot of the Airy beam is well integrated into the vortex part, forming a series of concentric rings similar to the Airy vortex beam, and its vortex state can be further confirmed in the phase field. The beam combines the singularity properties of vortices and the self-healing properties of Airy beams, which makes it promising for future laser applications.
Finally, under the control of adding HOCP to the Bessel vortex beam, we analyze the polygonal Airy vortex-like beam generated by introducing HOCP. Figure 9 shows the light field schematic diagrams of the xoz plane at the focus position, corresponding to the polygonal Airy vortex-like beams generated after the introduction of HOCP from the 3rd order to the 6th order, respectively. The generation of triangular, quadrilateral, pentagonal and hexagonal Airy vortex-like beams is well achieved at the focal position. It can be clearly analyzed in Figure 9 that the Airy beam has been well integrated into the polygonal Bessel beam, and the beam still retains the properties of the Airy beam. This design provides a reference value for the study of special vortex beams in the optical field.

Conclusions
In summary, we simulated an all-dielectric multifunctional metasurface with high transmittance by FDTD and realized three different functions of generating polygonal Bessel vortex beams, Airy beams and polygonal Airy vortex-like beams. The generation of polygons mainly depends on a new phase structure, namely, HOCP. We confirmed xoy plane phase field of the beam generated by RCP incidence, LP incidence and LCP incidence were simulated without HOCP modulation.
Finally, under the control of adding HOCP to the Bessel vortex beam, we analyze the polygonal Airy vortex-like beam generated by introducing HOCP. Figure 9 shows the light field schematic diagrams of the xoz plane at the focus position, corresponding to the polygonal Airy vortex-like beams generated after the introduction of HOCP from the 3rd order to the 6th order, respectively. The generation of triangular, quadrilateral, pentagonal and hexagonal Airy vortex-like beams is well achieved at the focal position. It can be clearly analyzed in Figure 9 that the Airy beam has been well integrated into the polygonal Bessel beam, and the beam still retains the properties of the Airy beam. This design provides a reference value for the study of special vortex beams in the optical field. Nanomaterials 2023, 13, x FOR PEER REVIEW 10 of 13 Figure 8. The (a1-a3) xoz plane light field, (b1-b3) xoy plane light field and (c1-c3) xoy plane phase field of the beam generated by RCP incidence, LP incidence and LCP incidence were simulated without HOCP modulation.
Finally, under the control of adding HOCP to the Bessel vortex beam, we analyze the polygonal Airy vortex-like beam generated by introducing HOCP. Figure 9 shows the light field schematic diagrams of the xoz plane at the focus position, corresponding to the polygonal Airy vortex-like beams generated after the introduction of HOCP from the 3rd order to the 6th order, respectively. The generation of triangular, quadrilateral, pentagonal and hexagonal Airy vortex-like beams is well achieved at the focal position. It can be clearly analyzed in Figure 9 that the Airy beam has been well integrated into the polygonal Bessel beam, and the beam still retains the properties of the Airy beam. This design provides a reference value for the study of special vortex beams in the optical field.

Conclusions
In summary, we simulated an all-dielectric multifunctional metasurface with high transmittance by FDTD and realized three different functions of generating polygonal Bessel vortex beams, Airy beams and polygonal Airy vortex-like beams. The generation of polygons mainly depends on a new phase structure, namely, HOCP. We confirmed

Conclusions
In summary, we simulated an all-dielectric multifunctional metasurface with high transmittance by FDTD and realized three different functions of generating polygonal Bessel vortex beams, Airy beams and polygonal Airy vortex-like beams. The generation of polygons mainly depends on a new phase structure, namely, HOCP. We confirmed from theoretical calculations and numerical simulations that the metasurface can well shape general vortex beams into arbitrary polygons by controlling the order of the HOCP. Additionally, it is further discovered that the HOCP only changes the light field intensity distribution but not the OAM spectrum. All generated beams have relatively narrow beam widths and extended undiffracted propagation distances. This simple and general approach provides another feasible way to realize the design of multifunctional components. Notably, the functions generated by the metasurface do not interfere with each other, which means the system possesses potential applications in the fields of particle trapping, 3D optical tweezers, parallel laser printing and multiplex fluorescence imaging.