Generation and Superposition of Perfect Vortex Beams in Terahertz Region via Single-Layer All-Dielectric Metasurface

Terahertz (THz) orbital angular momentum (OAM) technology provides promising applications in future wireless communication with large bandwidth and high capacity. However, the ring radius of the conventional THz vortex beam is related to the topological charge, limiting the co-propagation of multiple OAM modes in the THz communication systems. Although the perfect vortex beam (PVB) based on traditional methods can solve this problem, they are usually bulky and unstable. Here, we demonstrate two PVB generators based on a single all-dielectric metasurface to obtain polarization-independent PVB and spin multiplexed PVB, respectively. The former regulates the propagation phase by using isotropic unit cells; the latter simultaneously manipulates the propagation and geometric phase to achieve the spin-decoupled phase control by arranging anisotropic unit cells. In addition, we also demonstrate the stable generation of a perfect Poincaré beam with arbitrary polarization and phase distribution on a hybrid-order Poincaré Sphere via a spin-decoupled metasurface, which is achieved by the linear superposition of two PVBs with orthogonal circular polarizations. The proposed scheme provides a compact and efficient platform for the generation and superposition of PVBs in THz region, and will speed up the progress of THz communication systems, complex light field generation, and quantum information sciences.


Introduction
It is well known that light has both spin angular momentum (SAM) and orbital angular momentum (OAM), which are manifested as circular polarization and azimuthal phase of a light beam, respectively [1,2]. Owing to the annular intensity distribution and the particular helical phase structure, light with OAM (namely vortex beams) has been widely used in optical communication [3], optical trapping [4], metrology [5], and quantum memories [6]. Terahertz (THz) waves located between far-infrared and microwave bands with rich frequency resources and have not been utilized extensively [7]. The related wireless THz communication has been recognized as one of the potential key technologies for high-speed communication in the future due to its great performance in high rate, broad bandwidth, and strong anti-interference [8,9]. The combination of THz waves and vortex beams will further extend communication capacity because different OAM modes are strictly orthogonal, and the multiple OAM states of vortex beams can be used to transmit different data [10,11]. However, the ring radiuses of these conventional vortex beams strongly depend on their topological charges, so it is not easy to combine multiple OAM beams into a single fiber used for multiplexed communications simultaneously [12,13].
A similar situation occurred in the vector vortex beams (VVBs) domain. Such spiral phase distribution and spatial inhomogeneous polarization state provide more interaction between light and matter, and have been utilized in high-resolution lithography [14], example, it cannot be utilized to generate PVBs composed of arbitrarily circularly polarized vortex beams because there will always be one diverging circularly polarized output due to reversed phase sign. Therefore, to our best knowledge, a compact setup capable of generating an arbitrary PVB and further superposition to obtain PPBs in THz region is highly desirable and remains unexplored.
Here, two kinds of single-layer all-silicon metasurfaces have been designed to generate PVBs and PPBs in THz region. One generates polarization-independent PVBs by the propagation phase of cross-shaped unit cells. The other, made up of rectangular-shaped silicon pillars, can simultaneously employ the geometric and dynamic phase to generate spin-decoupled PVBs. Therefore, the radius and topological charges of the PVB can be conveniently varied by switching the incident spin states. Most importantly, two completely different PVBs can be mixed by a spin-decoupled metasurface. All orders of PPBs with different phase and polarization distributions can be generated by choosing the polarization states of incident light and designing topological charges of output PVBs. Thus, our work provides a potential compact platform for accelerating THz communication systems, complex light field generation, and quantum information sciences.

Polarization-Independent Metasurface
We design a metasurface that can generate polarization-independent PVBs. As shown in Figure 1a,b, the incident light with arbitrary polarization through the metasurface can produce a vortex beam such that radius is irrelevant to the topological charge. Theoretically, the complex amplitude of the PVB can be described as [24]: where (ρ, θ) represent the polar coordinates, δ(·) is the Dirac function in the polar coordinate system, l is the topological charge, and ρ 0 is the radius of bright annular intensity. Due to the Dirac function only existing in the ideal state, this perfect vortex model is not accessible in practice. According to the reported work, PVBs can be generated by taking the Fourier transform of an ideal higher-order Bessel beam [26], which can be presented in the following form, E B (r, z) = J l (k r ρ) exp(ilθ + ik z z) Here, J l is the first kind of l-th order Bessel function; k = K 2 r + K 2 z = 2π λ , (k r , k z ) are the radial and longitudinal wave vectors, respectively, and λ is the incident wavelength. However, it is not easy to produce ideal Bessel beams experimentally. In most practical cases, the accessible Bessel model is the Bessel-Gaussian (BG) beam. Hence, the PVB beam is typically generated from the Fourier transform of the BG beam, which is usually achieved by a cascading optical setup, including a spiral phase plate, an axicon lens, and a Fourier transform lens [44]. Specifically, the generation of PVBs has two steps. First, a spiral phase plate and an axicon are used to convert the Gaussian beam to the corresponding BG beam. Second, a Fourier lens is utilized to transform the BG beam into PVBs. Due to the unprecedented light manipulating ability, a single-layer metasurface can work as a spiral phase plate, an axicon, and a Fourier transformation lens simultaneously, as shown in Figure 1c. The total phase profile can be expressed as: ϕ axicon (x, y) = −2π where x and y represent the position coordinates of the unit cell in the metasurface concerning the center of the origin, Equation (4) represents the spiral phase profile of an optical vortex with the topological charge l. Equation (5) shows an axicon phase profile, and d is the axicon period relating to the ring radius of the PVB. Equation (6) where x and y represent the position coordinates of the unit cell in the metasurface concerning the center of the origin, Equation (4) represents the spiral phase profile of an optical vortex with the topological charge l. Equation (5) shows an axicon phase profile, and d is the axicon period relating to the ring radius of the PVB. Equation (6) introduces the phase profile of the Fourier transform lens with the focal length of f, and λ is the working wavelength. In order to realize the design requirement, we used the finite-difference time-domain (FDTD) method to optimize the geometric parameters of the unit cell. Single-layer plasmonic structures, depending on metallic materials with varying geometries, have been utilized to construct metasurfaces. However, they suffer from low efficiency due to the inherent ohmic loss and weak coupling between the incident and cross-polarized fields [45]. Although multilayer plasmonic structures can increase the efficiency, they also increase fabrication difficulties [46]. In addition, plasmonic structures are sensitive to the polarization of incident light, making them unsuitable for practical applications. It has been proved that the most straightforward approach to achieve polarization-independent performance is to find the unit cell with 4-fold symmetry [47]. Therefore, as shown in Figure 2a, the cross-shaped unit cell with the lattice constant of P = 130 μm, the height of H = 300 μm, and the exact width W = 50 μm are selected at the working frequency 0.75 THz. High-resistance silicon (ρ > 5000 Ω m) has been chosen as the designed material due to the negligible absorption and very low dissipation for terahertz wave as well as mature etching techniques, and its dielectric constant is set as ε = 11.9 [48,49].
By scanning the length (L) of the unit cell from 40 μm to 125 μm, both the transmittance and phase shift have been calculated, as shown in Figure 2b. Moreover, the In order to realize the design requirement, we used the finite-difference time-domain (FDTD) method to optimize the geometric parameters of the unit cell. Single-layer plasmonic structures, depending on metallic materials with varying geometries, have been utilized to construct metasurfaces. However, they suffer from low efficiency due to the inherent ohmic loss and weak coupling between the incident and cross-polarized fields [45]. Although multilayer plasmonic structures can increase the efficiency, they also increase fabrication difficulties [46]. In addition, plasmonic structures are sensitive to the polarization of incident light, making them unsuitable for practical applications. It has been proved that the most straightforward approach to achieve polarization-independent performance is to find the unit cell with 4-fold symmetry [47]. Therefore, as shown in Figure 2a, the cross-shaped unit cell with the lattice constant of P = 130 µm, the height of H = 300 µm, and the exact width W = 50 µm are selected at the working frequency 0.75 THz. High-resistance silicon (ρ > 5000 Ω m) has been chosen as the designed material due to the negligible absorption and very low dissipation for terahertz wave as well as mature etching techniques, and its dielectric constant is set as ε = 11.9 [48,49].
By scanning the length (L) of the unit cell from 40 µm to 125 µm, both the transmittance and phase shift have been calculated, as shown in Figure 2b. Moreover, the inset illustrates the side views of the simulated magnetic amplitude distributions when L is fixed as 120 µm. The magnetic fields are mainly localized inside the silicon pillars, indicating that the unit cell can be seen as a truncated low-loss waveguide that works as a weakly coupled low-quality factor Fabry-Perot resonator [50]. Next, as shown in Figure 2c, eight units are chosen to realize the complete phase control with high transmittance. Their phase shifts cover the 2π range, and the transmittances are maintained over 70%. Figure 2d shows phase wavefronts under the x-polarized wave incident at the pillar height H of 330 µm and the length L of 40, 82, and 120 µm, respectively. Clearly, the relative phase difference between these pillars at the output port (at 330 µm) corresponds to π, and 2π, respectively. a weakly coupled low-quality factor Fabry-Perot resonator [50]. Next, as shown in Figure  2c, eight units are chosen to realize the complete phase control with high transmittance. Their phase shifts cover the 2π range, and the transmittances are maintained over 70%. Figure 2d shows phase wavefronts under the x-polarized wave incident at the pillar height H of 330 μm and the length L of 40, 82, and 120 μm, respectively. Clearly, the relative phase difference between these pillars at the output port (at 330 μm) corresponds to π, and 2π, respectively.

Spin-Decoupled Metasurface
We design a spin-decoupled metasurface to generate two completely different PVBs only by varying the incident spin state. As shown in Figure 3a, left-circularly polarized (LCP) incident light through the metasurface can be transformed into a PVB with opposite-handedness, and the same metasurface can also transform right-circularly polarized (RCP) incident light into a different PVB with LCP state. Theoretically, the device must provide two independent spatial phase profiles corresponding to LCP and RCP, which cannot be achieved by using a single geometric phase or propagation phase.
For an anisotropic unit cell, the corresponding Jones matrix is: cos sin sin cos sin cos , sin cos sin cos cos sin where = ⅇ , = ⅇ are complex amplitude, and represent the transmission amplitude and the phase delay while the incident light is along with the ⅈ-

Spin-Decoupled Metasurface
We design a spin-decoupled metasurface to generate two completely different PVBs only by varying the incident spin state. As shown in Figure 3a, left-circularly polarized (LCP) incident light through the metasurface can be transformed into a PVB with oppositehandedness, and the same metasurface can also transform right-circularly polarized (RCP) incident light into a different PVB with LCP state. Theoretically, the device must provide two independent spatial phase profiles corresponding to LCP and RCP, which cannot be achieved by using a single geometric phase or propagation phase.
For an anisotropic unit cell, the corresponding Jones matrix is: where t x = T x e iϕ x , t y = T y e iϕ y are complex amplitude, T i and ϕ i represent the transmission amplitude and the phase delay while the incident light is along with the i-polarized wave with i (x, y), respectively. R(θ) = cos θ − sin θ sin θ cos θ is the rotation matrix with the rotation angle θ. When the metasurface is illuminated by circularly polarized light as E in = the output electric field can be expressed as: Here 1 i and 1 −i represent LCP and RCP, respectively. The transmitted light can be divided into two parts: the co-polarized component without any phase regulation and the cross-polarized component with conjugated phase modulation ±2θ. Combining the orientation-dependent geometric phase with the dimension-dependent propagation phase effectively will address the spin-locked limitation of geometry metasurfaces [51,52]. Consider the unit cell is lossless (T x = T y = 1) with π phase difference Nanomaterials 2022, 12, 3010 , indicating the pillars working as a half waveplate. The output electric field can then be expressed as: Assuming the expected phase distributions of the spin-decoupled are ϕ 1 (x, y) and ϕ 2 (x, y). From the Equation (9), we can derive that: where ϕ x (x, y) represents the propagation phase, which is influenced by the material and geometry parameters of the unit pillars; ±2θ(x, y) is the geometric phase, which is only determined by the rotation angle of the unit cell. Therefore, the expected phase distributions ϕ 1 (x, y) and ϕ 2 (x, y) can be satisfied simultaneously by designing the geometry and rotation angle of the unit structure, as shown in Figure 3b.
Here [ 1 ⅈ ] and [ 1 −ⅈ ] represent LCP and RCP, respectively. The transmitted light can be divided into two parts: the co-polarized component without any phase regulation and the cross-polarized component with conjugated phase modulation ±2 . Combining the orientation-dependent geometric phase with the dimension-dependent propagation phase effectively will address the spin-locked limitation of geometry metasurfaces [51,52]. Consider the unit cell is lossless ( = = 1 ) with π phase difference ( | ( , ) − ( , ) = π|), indicating the pillars working as a half waveplate. The output electric field can then be expressed as: Assuming the expected phase distributions of the spin-decoupled are 1 ( , ) and 2 ( , ). From the Equation (9), we can derive that: where ( , ) represents the propagation phase, which is influenced by the material and geometry parameters of the unit pillars; ±2 ( , ) is the geometric phase, which is only determined by the rotation angle of the unit cell. Therefore, the expected phase distributions 1 ( , ) and 2 ( , ) can be satisfied simultaneously by designing the geometry and rotation angle of the unit structure, as shown in Figure 3b. The designed total phase of the single-layer metasurface simultaneously contains the geometric phase and the dynamic phase, which are only related to geometry parameters and the rotation angle of the unit structure, respectively.
To construct the metasurface, the sub-wavelength rectangular-shaped silicon pillar has been chosen, as shown in Figure 4a. The simulation is based on the FDTD method to obtain the pillars' transmittance and phase delay database under the x-and y-polarized incident waves. The length L and width W varied from 35 μm to 135 μm with the fixed lattice constant of P = 140 μm and the height of H = 330 μm at the working frequency of 0.75 THz. Figure 4b,c show the transmittance as the function of the pillar sizes in L and To construct the metasurface, the sub-wavelength rectangular-shaped silicon pillar has been chosen, as shown in Figure 4a. The simulation is based on the FDTD method to obtain the pillars' transmittance and phase delay database under the x-and y-polarized incident waves. The length L and width W varied from 35 µm to 135 µm with the fixed lattice constant of P = 140 µm and the height of H = 330 µm at the working frequency of 0.75 THz. Figure 4b,c show the transmittance as the function of the pillar sizes in L and W, respectively. The corresponding simulated phase shifts (ϕ x , ϕ y ) are illustrated in Figure 4d,e. Both ϕ x and ϕ y cover 2π degrees in the scanning range, and all the transmittances are more than 70%. Therefore, by carefully selecting L and W from the database, a nearly arbitrary phase combination (ϕ x , ϕ y ) with high transmittance can be obtained for the silicon pillar. For our design, nine different silicon pillars are chosen, and the corresponding transmittances, phase shifts, and polarization conversion efficiencies (PCE, defined as the intensity ratio between the converted spin component and the total transmitted beam) are shown in Figure 4f. Phase shifts ϕ x (red square) and ϕ y (green square) are linearly increased with a step of 2π/9. The corresponding transmittance T x (red sphere) and T y (green sphere) are more than 70%. Moreover, the PCE (orange column) is maintained above 70%. Most importantly, the phase differences between ϕ x and ϕ y are close to π, and the transmittance T x and T y are almost equal, indicating these selected pillars are working well as the halfwave plates at 0.75 THz. The specific geometric parameters and PCE of the nine selected silicon pillars are listed in Table 1. transmittances are more than 70%. Therefore, by carefully selecting L and W from the database, a nearly arbitrary phase combination ( , ) with high transmittance can be obtained for the silicon pillar. For our design, nine different silicon pillars are chosen, and the corresponding transmittances, phase shifts, and polarization conversion efficiencies (PCE, defined as the intensity ratio between the converted spin component and the total transmitted beam) are shown in Figure 4f. Phase shifts (red square) and (green square) are linearly increased with a step of 2π/9. The corresponding transmittance (red sphere) and (green sphere) are more than 70%. Moreover, the PCE (orange column) is maintained above 70%. Most importantly, the phase differences between and are close to π, and the transmittance and are almost equal, indicating these selected pillars are working well as the half-wave plates at 0.75 THz. The specific geometric parameters and PCE of the nine selected silicon pillars are listed in Table 1.

Superpositions of PVBs by a Spin-Decoupled Metasurface
Under the paraxial approximation, the vectorial field of a monochromatic Poincaré beam can be decomposed into the superposition of two orthogonal circularly polarized vortex beams with different coefficients and as described [22]: where | ⟩ = ⅇ ( + ⅈ )/√2, | ⟩ = ⅇ ( − ⅈ )/√2 are the RCP and LCP vortex beams with topological charges of and m, respectively, while ex and ey are the unit vectors along the x and y axes.
is the azimuthal angle in the polar coordinates. In order

Superpositions of PVBs by a Spin-Decoupled Metasurface
Under the paraxial approximation, the vectorial field of a monochromatic Poincaré beam can be decomposed into the superposition of two orthogonal circularly polarized vortex beams with different coefficients ψ l N and ψ m s as described [22]: where |R l = e il ϕ (e x + ie y )/ √ 2 , |L m = e imϕ (e x − ie y / √ 2 are the RCP and LCP vortex beams with topological charges of l and m, respectively, while e x and e y are the unit vectors along the x and y axes. ϕ is the azimuthal angle in the polar coordinates. In order to obtain PPBs that radius is irrelevant to topological charges, Equation (11) can be further deduced as: Here, |POV R , l = e iϕ total e x + ie y / √ 2and|POV L , m = e iϕ total e x − ie y / √ 2 represents the RCP and LCP PVBs with topological charges of l and m, ϕ total represent the phase profile of PVBs obtained from Equation (3).
According to Equation (12), two essential parts, orthogonal circularly polarized PVBs and coefficients ψ l N and ψ m s , should be resolved for generating arbitrary PPBs. In order to describe the process more clearly, we use the PS to describe scalar fields of incident light Nanomaterials 2022, 12, 3010 8 of 16 and the HyOPS to describe the vectorial field of output PPBs. As shown in Figure 5a, all the possible states of fundamental polarization can be mapped onto the PS. Because the polarization state of the scalar field can be represented with orthogonal circular polarization, the coefficients ψ l N and ψ m s vary with the incident polarization. As shown in Figure 5b, the spin-decoupled metasurface can generate two completely different PVBs with orthogonal circularly polarization. If the axicon period and the focal length are identical, arbitrary PPBs can be generated by superposing the two circularly polarized PVBs on a HyOPS. As shown in Figure 5c, each point on the surface of HyOPS with spherical coordinates (α, β) can indicate a PPB, where α ∈ (0, 2π), β ∈ (0, π). The two poles of the HyOPS represent two PVBs with opposite circular polarization and topological charges of l and m, respectively. Points on the equator of HyOPS represent a superposition of equal intensities and different phase differences of the opposite circular polarization states. For example, the points (0, π/2), (π, π/2) represent the radially and azimuthally vector beam, respectively.
According to Equation (12), two essential parts, orthogonal circularly polarized PVBs and coefficients and , should be resolved for generating arbitrary PPBs. In order to describe the process more clearly, we use the PS to describe scalar fields of incident light and the HyOPS to describe the vectorial field of output PPBs. As shown in Figure 5a, all the possible states of fundamental polarization can be mapped onto the PS. Because the polarization state of the scalar field can be represented with orthogonal circular polarization, the coefficients and vary with the incident polarization. As shown in Figure 5b, the spin-decoupled metasurface can generate two completely different PVBs with orthogonal circularly polarization. If the axicon period and the focal length are identical, arbitrary PPBs can be generated by superposing the two circularly polarized PVBs on a HyOPS. As shown in Figure 5c, each point on the surface of HyOPS with spherical coordinates ( , ) can indicate a PPB, where ∈ (0,2π), ∈ (0, π). The two poles of the HyOPS represent two PVBs with opposite circular polarization and topological charges of and , respectively. Points on the equator of HyOPS represent a superposition of equal intensities and different phase differences of the opposite circular polarization states. For example, the points (0, π/2), (π, π/2) represent the radially and azimuthally vector beam, respectively. For the PPBs on a HyOPS, the polarization state can be mapped by representing the Stokes parameters in the spherical Cartesian coordinates, which are defined as follows [20]: For the PPBs on a HyOPS, the polarization state can be mapped by representing the Stokes parameters in the spherical Cartesian coordinates, which are defined as follows [20]: Here Φ = arg ψ l N ) − arg ψ l N is the phase difference between the two PVBs. ψ l N 2 and ψ m S 2 are the intensity of RCP and LCP PVB, respectively. The concepts of polarization order P = (l − m)/2 and the topological Pancharatnam charge l P = (l + m)/2 are also used to characterize PPBs [22]. This can be explained by deducting Equation (12) as: Equation (17) indicates that the polarization state of PPBs is only determined by the parameters within the brackets, that is, the polarization order P and the HyOPS coordinate(α, β). The topological Pancharatnam charge l P characterizes the phase distribution of the PPBs, which can be represented as e iϕl P .

Results and Discussion
Based on the above design principles, we designed and numerically simulated several metasurface devices (diameter is 14 mm) to demonstrate the practicability of our proposed design by employing the FDTD method. The perfectly matching layers (PMLs) were used in the x, y, and z directions. Plane-wave sources were used in all simulations.

Polarization-Independent PVBs Generator
The generation of polarization-independent PVBs has been investigated. The design phase distribution has been described in Equation (3). Figure 6a displays the intensity profiles in the x-z plane of the generated PVB with the parameters of d = 4 mm and l = 2 at 0.75 THz. Figure 6b plots the corresponding cross-section images of the generated PVB at different longitudinal positions. The THz beam passing through the metasurface presents an annular ring profile near the focal plane (z = 14 mm). Moreover, the intensity profile of the ring maintains a circular shape, but the radius increases with increasing propagation distance.
| | are the intensity of RCP and LCP PVB, respectively. The concepts of polarization order = ( − )/2 and the topological Pancharatnam charge = ( + )/2 are also used to characterize PPBs [22]. This can be explained by deducting Equation (12) Equation (17) indicates that the polarization state of PPBs is only determined by the parameters within the brackets, that is, the polarization order and the HyOPS coordinate ( , ) . The topological Pancharatnam charge characterizes the phase distribution of the PPBs, which can be represented as ⅇ ⅈ .

Results and Discussion
Based on the above design principles, we designed and numerically simulated several metasurface devices (diameter is 14 mm) to demonstrate the practicability of our proposed design by employing the FDTD method. The perfectly matching layers (PMLs) were used in the x, y, and z directions. Plane-wave sources were used in all simulations.

Polarization-Independent PVBs Generator
The generation of polarization-independent PVBs has been investigated. The design phase distribution has been described in Equation (3). Figure 6a displays the intensity profiles in the x-z plane of the generated PVB with the parameters of d = 4 mm and l = 2 at 0.75 THz. Figure 6b plots the corresponding cross-section images of the generated PVB at different longitudinal positions. The THz beam passing through the metasurface presents an annular ring profile near the focal plane (z = 14 mm). Moreover, the intensity profile of the ring maintains a circular shape, but the radius increases with increasing propagation distance.   Figure 7b shows the corresponding cross-sections of the normalized intensity profiles extracted along the x-direction (solid line) and y-direction (dashed line). The consistent cross-sectional intensity profiles show that the generated vortex beams are insensitive to the topological charge, and the phase distribution is consistent with the designed topological charge. The full-wave simulation generation   Figure 7b shows the corresponding cross-sections of the normalized intensity profiles extracted along the x-direction (solid line) and y-direction (dashed line). The consistent cross-sectional intensity profiles show that the generated vortex beams are insensitive to the topological charge, and the phase distribution is consistent with the designed topological charge. The full-wave simulation generation efficiencies (defined as the ratio of the PVB's intensity to the intensity of the incident beam) of all polarization-independent metasurfaces are more than 49%. The polarization-independent properties and the relationship between the ring radius and the axicon period are analyzed in Sections S1 and S2 of the Supplemental Materials. Section S3 of the Supplemental Materials characterizes the working bandwidth of the polarization-independent PVBs generator. efficiencies (defined as the ratio of the PVB's intensity to the intensity of the incident beam) of all polarization-independent metasurfaces are more than 49%. The polarizationindependent properties and the relationship between the ring radius and the axicon period are analyzed in Sections S1 and S2 of the Supplemental Materials. Section S3 of the Supplemental Materials characterizes the working bandwidth of the polarizationindependent PVBs generator.

Spin-Decoupled PVBs Generator
According to theoretical analysis of a spin-decoupled PVBs generator, a single metasurface can generate two distinct PVBs. The corresponding phase distribution ( 1 , 2 ) can be expressed as follows: x y x y l y d f where λ is the incident wavelength, li, fi, and di represent corresponding parameters for the incident RCP (i = 1) and LCP (i = 2) light, respectively. Two metasurface devices are designed based on Equations (10) and (18) to demonstrate the above function. The first metasurface is designed with l1 = 2, l2 = 3, d1 = d2 = 4 mm, and f1 = f2 = 14 mm at 0.75 THz. As shown in Figure 8a,c, when the incident light is converted from RCP to LCP, the transmitted light is converted into a corresponding cross-polarized component, and the ring radius does not change significantly, but the topological charge l changes from 2 to 3, which is consistent with our design. In addition, we can also control the ring radius of generated PVBs by switching the incident spin state. The second device with the parameters l1 = l2 = 3, d1 = 4 mm, d2 = 6 mm, and f1 = f2 =14 mm at 0.75 THz is also designed. As shown in Figure 8b,d, the topological charge value l does not change, but the ring radius is inversely proportional to the axicon periods d when the incident light is

Spin-Decoupled PVBs Generator
According to theoretical analysis of a spin-decoupled PVBs generator, a single metasurface can generate two distinct PVBs. The corresponding phase distribution (ϕ 1 , ϕ 2 ) can be expressed as follows: where λ is the incident wavelength, l i , f i , and d i represent corresponding parameters for the incident RCP (i = 1) and LCP (i = 2) light, respectively. Two metasurface devices are designed based on Equations (10) and (18) to demonstrate the above function. The first metasurface is designed with l 1 = 2, l 2 = 3, d 1 = d 2 = 4 mm, and f 1 = f 2 = 14 mm at 0.75 THz. As shown in Figure 8a,c, when the incident light is converted from RCP to LCP, the transmitted light is converted into a corresponding cross-polarized component, and the ring radius does not change significantly, but the topological charge l changes from 2 to 3, which is consistent with our design. In addition, we can also control the ring radius of generated PVBs by switching the incident spin state. The second device with the parameters l 1 = l 2 = 3, d 1 = 4 mm, d 2 = 6 mm, and f 1 = f 2 =14 mm at 0.75 THz is also designed. As shown in Figure 8b,d, the topological charge value l does not change, but the ring radius is inversely proportional to the axicon periods d when the incident light is converted from RCP to LCP. It is worth noting that the transmitted light still carries with a co-polarized component, as shown in Figure 8a,b, because the selected unit cells are not a perfect half-wave plate. The full-wave simulation generation efficiencies of the spindecoupled metasurfaces were more than 53%, which can be improved by selecting materials with higher transmittance and a more careful phase scanning operation with smaller step size (more details on working bandwidth of the meta-device seen also in Section S4 of the Supplemental Materials). Therefore, the spin-decoupled and polarization-independent PVB meta-devices will pave the way for developing a miniaturized and convenient platform for OAM-related applications in the terahertz band, especially in improving the transmission capacity of terahertz communications.
decoupled metasurfaces were more than 53%, which can be improved by selecting materials with higher transmittance and a more careful phase scanning operation with smaller step size (more details on working bandwidth of the meta-device seen also in Section S4 of the Supplemental Materials). Therefore, the spin-decoupled and polarization-independent PVB meta-devices will pave the way for developing a miniaturized and convenient platform for OAM-related applications in the terahertz band, especially in improving the transmission capacity of terahertz communications.

Superpositions of Perfect Vortex Beams
According to the theoretical analysis of generating PPBs in Section 2.3, the topological charges (l, m) of orthogonal circularly polarized PVBs and the corresponding coefficients and can determine a PPB. Specifically, different topological charges (l, m) satisfied by spin-decoupled metasurface will determine a unique HyOPS, different coefficients and will define the position of a PPB on this HyOPS by calculating Stokes parameters. In order to verify the above principle, three spin-decoupled metasurfaces with the expected parameters d = 4 mm, and f =14 mm at 0.75 THz have been designed. The relevant polarization order and topological Pancharatnam charge are 1 = ( 1 − 1 )/2 = −1 , 1 = ( 1 + 1 )/2 =0, 2 = ( 2 − 2 )/2 = −2 , 2 = ( 2 + 2 )/2 = 0, and 3 = ( 3 − 3 )/2 = −1 , 3 = ( 3 + 3 )/2 = 2. Five points with their coordinates on the HyOPS are selected, given by the red dots in Figure 9a, and the corresponding incident polarization states are represented by the red points located on the PS shown in Figure 9b. Figure 9c shows the intensity patterns of three designed metasurfaces in the x-y plane at the focus point. The annular intensity patterns are obtained by a horizontal linear polarizer to analyze the polarization property. Because the

Superpositions of Perfect Vortex Beams
According to the theoretical analysis of generating PPBs in Section 2.3, the topological charges (l, m) of orthogonal circularly polarized PVBs and the corresponding coefficients ψ l N and ψ m s can determine a PPB. Specifically, different topological charges (l, m) satisfied by spin-decoupled metasurface will determine a unique HyOPS, different coefficients ψ l N and ψ m s will define the position of a PPB on this HyOPS by calculating Stokes parameters. In order to verify the above principle, three spin-decoupled metasurfaces with the expected parameters d = 4 mm, and f = 14 mm at 0.75 THz have been designed. The relevant polarization order P and topological Pancharatnam charge l P are P 1 = (l 1 − m 1 )/2 = −1, l P1 = (l 1 + m 1 )/2=0, P 2 = (l 2 − m 2 )/2 = −2, l P2 = (l 2 + m 2 )/2 = 0, and P 3 = (l 3 − m 3 )/2 = −1, l P3 = (l 3 + m 3 )/2 = 2. Five points with their coordinates on the HyOPS are selected, given by the red dots in Figure 9a, and the corresponding incident polarization states are represented by the red points located on the PS shown in Figure 9b. Figure 9c shows the intensity patterns of three designed metasurfaces in the x-y plane at the focus point. The annular intensity patterns are obtained by a horizontal linear polarizer to analyze the polarization property. Because the polarization order P determines the number of polarization rotations per round trip, the transmitted intensity patterns through a horizontal linear polarizer will have 2|P| lobes.
The case with the same topological Pancharatnam charge l P and different polarization order P is firstly discussed. As shown in the first row and the second row in Figure 9c, the transmitted intensity patterns are split into two and four lobes, respectively, which satisfy the theoretical predictions of P 1 = −1, P 2 = −2. Then, the influence of the topological Pancharatnam charge l P is considered. The polarization order is fixed as P = −1, and varied with the topological Pancharatnam charge l P = 0, 2, respectively. The simulated annular intensity distributions in the first and third row show the same annular patterns, demonstrating the independence of the polarization state of PPBs from the topological Pancharatnam charge, matching well with the theoretical prediction of Equation (17). It should be noticed that the third row has a slight rotation of the measured intensity induced by the Gouy phase [53]. In principle, the Gouy phase is dependent on the topological charge, and the rotation of the intensity distribution around the beam axis occurs only when the topological charges |l| = |m|. Therefore, this approach provides a potential solution to measure the Gouy phase.
which satisfy the theoretical predictions of 1 = −1, 2 = −2. Then, the influence of the topological Pancharatnam charge is considered. The polarization order is fixed as = −1, and varied with the topological Pancharatnam charge = 0, 2, respectively. The simulated annular intensity distributions in the first and third row show the same annular patterns, demonstrating the independence of the polarization state of PPBs from the topological Pancharatnam charge, matching well with the theoretical prediction of Equation (17). It should be noticed that the third row has a slight rotation of the measured intensity induced by the Gouy phase [53]. In principle, the Gouy phase is dependent on the topological charge, and the rotation of the intensity distribution around the beam axis occurs only when the topological charges | | ≠ | |. Therefore, this approach provides a potential solution to measure the Gouy phase. To further characterize the capability of arbitrary PPBs generation, another five incident lights represented by the blue points on the Poincaré sphere are chosen, as shown in Figure 9b, where the two circular eigenstates have a phase difference of π/2. The corresponding PPBs with their coordinates on the HyOPS are expressed by the blue dots in Figure 9a, and the measured intensity patterns through a horizontal linear polarizer are shown in Figure 9d. Unsurprisingly, the introduced phase difference of π/2 between the To further characterize the capability of arbitrary PPBs generation, another five incident lights represented by the blue points on the Poincaré sphere are chosen, as shown in Figure 9b, where the two circular eigenstates have a phase difference of π/2. The corresponding PPBs with their coordinates on the HyOPS are expressed by the blue dots in Figure 9a, and the measured intensity patterns through a horizontal linear polarizer are shown in Figure 9d. Unsurprisingly, the introduced phase difference of π/2 between the two circular eigenstates results in the change of coefficients ψ l N and ψ m s , further causing the rotation of the intensity profiles and demonstrating the feasibility of PPBs by changing the polarization of incident light.
In addition, the Stokes parameters have also been used to verify inhomogeneous polarization states of the PPBs, which are given by [54]: where I 0 , I 90 , I 45 Figure 10a shows the calculated results of normalized Stokes parameters by choosing point III in Figure 9a for three kinds of PPBs. The first Stokes parameter S 0 , representing the light intensity, shows the same ring pattern consistent with our theoretical analysis. By analyzing Stokes parameters, S 1 and S 2 , the polarization orientation angle σ can be determined by: The corresponding polarization orientations (orange double arrows) and distributions of the PPBs are shown in Figure 10b. It can be seen that the polarization orientations rotate 2π and 4π per round trip for cases of P 1 = P 3 = −1, P 2 = −2, respectively, which is consistent to the definition of polarization order. It should be noticed that the spherical coordinate of point III is (0, π/2), which means the Stokes parameter S 3 is zero theoretically. The error appearing in Figure 10a mainly comes from two reasons: the phase difference of the metasurface between the actual phase and the theoretical phase, and the act that the selected unit cells are not a perfect half-wave plate. The full-wave simulation generation efficiencies of the PPBs (defined as the ratio of the PPB's intensity to the intensity of the incident beam) were more than 51%. These results and the analysis demonstrate that spin-decoupled metasurfaces can efficiently generate arbitrary PPBs on the HyOPS in THz region. Therefore, the terahertz PVB meta-devices designed here will also provide a potential compact platform for vector beam related applications in the terahertz band, such as reconstructing incident photons' frequency and polarization state [55].

Conclusions
In conclusion, we propose a general method to generate PVBs using a single-layer dielectric metasurface in THz regime. This is achieved by integrating the function of spiral

Conclusions
In conclusion, we propose a general method to generate PVBs using a single-layer dielectric metasurface in THz regime. This is achieved by integrating the function of spiral phase plate, axicon, and Fourier transformation lens into a single-layer metasurface. The ring radiuses of these generated vortex beams are independent of the topological charges and inversely proportional to the axicon periods. Meanwhile, the unit-cell cross-section geometries with 4-fold symmetry can make the generated PVBs independent of incident polarization. Moreover, to realize more freedom, a spin-decoupled metasurface is designed by simultaneously manipulating the dynamic phase and geometric phase of the unit cells. The metasurface can generate spin-decoupled PVBs with arbitrary radius or topological charges by changing the incident spin states. In addition, the generation of an arbitrary perfect Poincaré beam at any position on the surface of the HyOPS is also demonstrated by superposing two perfect vortex beams. Various perfect Poincaré beams were generated, and the results agreed well with the theory. This work will provide more opportunities and possibilities for THz communication, complex light field generation, and quantum information sciences.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/nano12173010/s1, Figure S1: (a) The intensity of the PVBs at the focus point under different polarized incident light (depicted by the red double arrow). (b) Corresponding normalized cross-section of the annular intensity distribution along x-direction (solid line) and y-direction (dashed line). (c) The relationship between the ring radius and the axicon period with the unified parameters l = 2, f = 14 mm. Figure S2. Characterization of the bandwidth of the polarization-independent PVBs generator. (a) The intensity of the PVBs at the focus point with a frequency ranging from 0.65 THz to 0.9 THz. (b-g) Corresponding normalized cross-section of the annular intensity distribution along x-direction (blue line) and y-direction (red line). Figure S3. Analysis of the bandwidth of the spin-decoupled PVBs generator. (a) The intensity of the PVBs at the focus point with a frequency ranging from 0.5 THz to 1.0 THz. (b-g) Corresponding normalized cross-section of the annular intensity distribution along x-direction (blue line) and y-direction (red line).