1. Introduction

Terahertz waves, characterized by low photon energy, strong penetration, and high spectral resolution, show broad prospects in fields such as wireless communication, biomedical imaging, and security screening [1,2]. Traditional methods for generating vortex beams suffer from issues like large device volume, wavelength limitations, and low integration [3,4]. As a two-dimensional artificial electromagnetic structure, metasurfaces can flexibly control phase, amplitude, and polarization through subwavelength unit cells [5-7], providing a new pathway for the efficient generation and control of vortex beams [8]. In recent years, metasurfaces have demonstrated great potential in achieving spin-to-orbital angular momentum conversion [9], structured light manipulation [10], and even quantum entanglement [11], opening new directions for the development of integrated photonic devices.

This paper designs a transmissive metasurface based on polyimide dielectric and copper patches, operating at a frequency of 0.1 THz. The metasurface employs a double-L-shaped unit structure, utilizes the geometric phase principle to convert circularly polarized light into vortex beams, and achieves multi-channel beam splitting and topological charge manipulation based on far-field scattering and coding metasurface theory. Simulation results show that this structure maintains high transmission efficiency across multiple frequency bands, including 0.1933-0.1982 THz and 0.2354-0.2386 THz, supporting the generation of multi-band, multi-channel vortex beams, thereby providing a design foundation for integrated terahertz photonic devices.

2. Fundamental Theoretical Principles

2.1 Generalized Snell's Law

The traditional law of refraction describes the propagation behavior of light at the interface between homogeneous media. When an artificially designed phase discontinuity is introduced at the interface, anomalous reflection and refraction phenomena occur. This principle can be described by the Generalized Snell's Law [12]. Metasurfaces achieve arbitrary wavefront manipulation by designing the phase response of unit cells at subwavelength scales.

As shown in Figure 1, assuming light travels from point A to point B, the paths can be divided into ADB and AEB. The optical phase at point D is Φ, and at point E is Φ+dΦ. According to Fermat's principle, the optical path lengths of these two paths should be equal, hence their phase difference should also be equal. Considering the existence of a phase gradient dΦ/dx at the interface, the Generalized Snell's Law for the transmission direction can be derived:

(1)

where θ_i and θ_t are the incident and transmission angles, respectively; n_i and n_t are the refractive indices of the media;λ is the wavelength. This paper primarily utilizes transmissive metasurfaces, controlling the direction of the outgoing beam by introducing continuous phase variation (phase gradient) between unit cells. If the metasurface achieves phase coverage over a length Γ, the phase gradient can be expressed as dΦ/dx=2π/Γ, simplifying the transmission angle to:

(2)

2.2 Pancharatnam-Berry Geometric Phase Principle

Pancharatnam and Berry independently pointed out that when circularly polarized light passes through a structure with anisotropic scattering characteristics and this structure rotates in-plane, the cross-polarized component of the outgoing light acquires an additional phase related to the rotation angle, known as the geometric phase.

Taking vertically incident left-handed circularly polarized (LCP) light as an example, if the unit cell possesses anisotropic transmission characteristics and satisfies certain amplitude and phase conditions, the transmitted right-handed circularly polarized (RCP) component will carry a phase factor ei2α, where α is the rotation angle of the unit cell. This means that by simply rotating the anisotropic unit, phase modulation from 0 to 2π can be introduced in the cross-polarized component, and this phase is wavelength-independent, depending solely on the geometric rotation angle. This principle underlies the phase manipulation achieved in this paper using rotating double-L-shaped structures [13,14].

 field enhances the interaction between biomolecules and the surface, thereby increasing the sensitivity and response efficiency of the metasurface in biosensing applications. Fig. 2 c) and Fig. 2 (d) illustrate the magnetic field energy distributions in the x-y and x-z planes at 0.7376 THz, respectively. These figures demonstrate that the magnetic field of the metasurface is primarily concentrated at the bottom of the elliptical nanocolumns and on the substrate between the elliptical structures. The distribution of the magnetic field differs from that of the electric field, and the areas of concentration for both fields are closely related to the geometric shape of the elliptical columns.

2.3 Far-Field Scattering and Coding Metasurfaces

Coding metasurfaces represent unit cells with different phase responses using binary digits arranged in specific sequences, enabling control over far-field scattering patterns [15]. Its principle is analogous to phased array antennas. For a 1-bit coding metasurface, it consists of two unit cells with a phase difference of exactly π, assigned binary codes "0" and "1", respectively. Similarly, a 2-bit coding metasurface consists of 4 unit cells with phase differences of π/2each, corresponding to binary codes "00", "01", "10", and "11". For a 3-bit metasurface, it comprises 8 types of unit cells, each with a phase difference of π/4. This digital design concept greatly simplifies the design and optimization process of metasurfaces and has led to new directions such as multifunctional complex-amplitude manipulation [16].

As shown in Figure 2, for a metasurface composed of m×n unit cells, when a plane wave is normally incident, its far-field scattering pattern F(θ,φ) can be expressed as the product of the unit cell structure factor and the array factor. By designing the coding sequence, i.e., the unit phase distribution φ(m,n), the direction, number, and shape of the scattered beams can be controlled. In particular, by designing coding periods Γx and Γy along orthogonal directions (x and y), beam deflection in three-dimensional space can be achieved, with the main lobe direction satisfying:

(3)

(4)

This principle provides the theoretical basis for the subsequent realization of multi-channel, multi-beam vortex beams.

1.  2.4 Phase Arrangement Principle for Vortex Beams

Linearly polarized light can be decomposed into two circularly polarized components with opposite handedness, each with a unique output direction. By designing the output directions for LCP and RCP, beam splitting in different directions can be achieved. Adding orbital angular momentum with different topological charges to the split beams in different directions enables multi-channel terahertz vortex beams.

When a plane wave is incident at a certain angle, the phase distribution ϕ(x,y) of the metasurface needs to satisfy:

(5)

Where the operating wavelength is λ, the azimuth angle is φ_i, and the elevation angle is θ_i. The phase distribution of each corresponding unit structure (x, y) is defined according to the vortex phase distribution, which needs to satisfy:where the operating wavelength is λ, the azimuth angle is φ_i, the elevation angle is θ_i. According to the definition of vortex phase distribution, the phase arrangement for each unit cell x,y needs to satisfy:

(6)

Where the working wavelength is λ, the azimuth angle is φi, the elevation angle is θi, and the corresponding phase distribution of each unit structure (x, y) is based on the vortex phase. The topological charge number of the vortex beam is represented by lij, jl represents RL or LR, and the vortex phase value at each point x, y is represented by ϕ_ijvortex(x,y). By superimposing the above two equations, the directional phase equation and the vortex phase equation, the phase distribution of the metasurface can be utilized to generate a multi-channel vortex beam distribution definition, which needs to satisf：

where l_ij represents the topological charge of the vortex beam, ij denotes RL or LR, and ϕ_ijvortex(x,y) represents the vortex phase value at each point x,y. By superposing the directional control phase equation and the vortex phase equation above, the phase distribution of the metasurface can be utilized to generate multi-channel vortex beams:

(7)

Since it is necessary to consider that an incident light can be decomposed into a left-handed circularly polarized light (LCP) and a right-handed circularly polarized light (RCP), and in order to achieve the outgoing light to be directed towards (θ₁, φ₁) and (θ₂, φ₂) respectively under the condition of incident linearly polarized light, adjustments need to be made to the phase distribution of the metasurface:Considering that an incident light can be decomposed into LCP and RCP components, and to achieve outgoing light directed towards θ₁, φ₁ and θ₂, φ₂ under linearly polarized incidence, the phase distribution of the metasurface needs to be adjusted:

(8)

3 Terahertz Metasurface Unit Cell Design

To achieve high-efficiency terahertz wavefront manipulation, this section designs and optimizes a double-L-shaped transmissive metasurface unit based on the geometric phase principle. With the operating frequency range fixed at 0.1 THz, it is necessary to determine the optimal unit period and the five parameters of the L-shaped structure. First, by parameter scanning of a single L-shape, the optimal parameter ratio for the long and short axes corresponding to the maximum transmission amplitude is determined.

As shown in Figure 3(a), the metasurface unit features an L-shaped structure. The unit uses low-loss, highly flexible polyimide as the dielectric substrate, with a thickness of h = 70 μm. The surface metal layer is copper with a thickness of 0.2 μm. The period of the L-shaped unit is P = 1449 μm, with long axis l₁ = 1071 μm, short axis l₂ = 430.5 μm, offset l₃ = 200.65 μm, l₄ = 512.4 μm, and line width l₅ = 86.1 μm Based on this, a unit period conforming to this ratio is selected, and the entire unit period and structure are scaled proportionally to observe the shifting trend of the maximum transmission amplitude peak, aiming to position it near 0.1 THz. Since the period significantly influences the unit's transmission amplitude, while the L-shaped metal film structure has a lesser effect, the period is prioritized in parameter selection. Using LCP incident light, Figure 3(b) clearly shows that as the unit period increases from P = 1035 μm to P = 1449 μm, the peak of maximum amplitude gradually shifts towards 0.1 THz. This roughly determines P = 1449 μm as the optimal period, where the transmission amplitude at 0.1 THz is maximum, approximately 0.4926.

4 Realization of Vortex Beams Based on Interleaved Arrangement

To verify whether the missing phase information in an interleaved pattern affects the metasurface's ability to generate vortex beams under the same array size, we designed two types of metasurfaces: one with coordinates both odd or both even, and another with coordinates not both odd or both even. Through these designs, we can compare the differences between the vortex beams generated by the two metasurfaces, aiding in a better understanding of the metasurface's working principle.

4.1 Vortex Beam Generation with Same Parity Coordinates

Based on the design principles above and using the L-shaped unit from the previous section, a metasurface with same parity coordinates was designed. As shown in Figure 4, a 30×30 same-parity metasurface is presented. This metasurface consists of 30×30 unit cells. When the coordinates of the unit cells are both odd or both even (e.g., [1,1]; [2,2]), the deflection angle is set to 0°, and the orientation angle is adjusted to 45°. The coordinates of other points are obtained and arranged according to the theoretical calculations mentioned earlier.

Under the condition of a center frequency of 0.1 THz, with LCP incident light along the negative z-axis direction onto this same-parity coordinate metasurface, the following results are obtained.

From Figure 5, it can be seen that the beam exhibits a spiral shape, and the phase intensity at the center of the main lobe is very low. Observing the phase distribution reveals that one full rotation achieves 2πphase coverage, generating a vortex beam with topological charge l=1.

4.2 Vortex Beam Generation with Different Parity Coordinates

Similarly, using the same principles and unit structure, a metasurface with different parity coordinates was designed, as shown in Figure 6, creating a 30×30 different-parity metasurface. This metasurface also consists of 30×30 unit cells. When the coordinates of two units are not both odd or both even (e.g., [1,2], [2,3]), their deflection angle is set to 0°, and the orientation angle is adjusted to 135°. The coordinates of other points are determined through the previous theoretical calculations and arranged accordingly.

For the different-parity metasurface, also at 0.1 THz frequency with LCP incident along the negative z-axis, the following results are obtained.

As shown in Figure 7, due to the superposition of two vortex beams, it is difficult to distinguish them clearly from the amplitude distribution, only roughly observing a vortex formed at the center. However, by examining the phase distribution, it can be seen that a vortex beam with topological charge l=1 is formed at the array center, and another vortex beam with topological charge l=2 exists slightly above the center.

Thus, we have proven that under the same array size, the missing phase information in an interleaved pattern does not prevent the metasurface from generating vortex beams. Furthermore, after superposition, vortex beams with different topological charges can still be generated without mutual interference. This verifies that a single metasurface can integrate the functionalities of multiple metasurfaces, generating vortex beams with different topological charges, thereby simplifying the fabrication process.

4.3 Superposition of Interleaved Metasurfaces

Since the two metasurfaces are interleaved, according to the vortex beam arrangement principle, we superimposed the same-parity metasurface and the different-parity metasurface to obtain a new composite metasurface of size 30×30.

The simulation results under LCP incidence at 0.1 THz are shown in Figure 9. Similar to before, due to superposition, the amplitude distribution is difficult to resolve clearly, but the phase distribution reveals a vortex beam with l=1 at the array center and another with l=2 slightly above the center.

This further confirms that interleaved missing phase information does not hinder vortex beam generation, and superposition can yield different topological charges without interference, validating the concept of integrating multiple functions into a single metasurface to simplify processes.

 5 Realization of Multi-Channel Vortex Beams

To manipulate vortex beams based on interleaved metasurfaces, the phase manipulation principle and vortex phase distribution can be employed by adjusting the deflection angles of unit cells to achieve multi-channel vortex beam generation. Since linearly polarized light can be decomposed into LCP and RCP components, independent wavefront manipulation for LCP and RCP can be designed separately, causing incident linearly polarized light to exit in different directions, forming beam splitting. Building on this, different topological charge OAM is further introduced into each output channel, ultimately realizing the generation and control of multi-channel terahertz vortex beams.

Based on the geometric phase manipulation mechanism, MATLAB programming was used to calculate the required rotation angles for each unit cell. According to formula (8), setting the elevation angle to 45°and azimuth angles to 135°and 45°, a 30×30 metasurface phase distribution matrix was generated. This matrix was imported via CST's COM interface for automatic modeling, with the dielectric layer as polyimide and the metal layer as copper. Figure 10(a) illustrates the designed metasurface phase distribution principle, while Figures 10(b)-(c) show the corresponding unit arrangement and array structure.

5.1 Left-Handed Circularly Polarized Light Incidence

In the CST time-domain solver, an LCP plane wave with a center frequency of 0.1 THz was used for normal incidence, with light incident along the negative z-axis direction (from the backside of the metal film) onto the metasurface array.

Under illumination by LCP (Note: text says RCP but context implies LCP), the metasurface generates four vortex beams with different topological charges. From the 3D amplitude distribution of the metasurface shown in Figure 11(a), four beams at distinct angles can be clearly distinguished. Figure 11(b) shows that at an elevation of 45° and azimuth of 135°, the main lobe exhibits a distinct spiral shape with very low phase intensity at the rotation center, generating a vortex beam with topological charge +1, achieving 360° phase coverage. Similarly, Figure 11(c) shows that at an elevation of 45° and azimuth of 45°, the main lobe clearly displays a spiral shape with very weak phase intensity at the center, generating a vortex beam with topological charge +2, achieving 720° phase coverage. The black circles in the 2D phase-amplitude distribution of the metasurface in Figure 11(d) indicate the generation of vortex beams with topological charges +2 and +1 at the respective angles.

This result demonstrates the feasibility of our theoretical design. Specifically, our design successfully controls the propagation and radiation behavior of electromagnetic waves, achieving efficient beam manipulation. The simulation results show the expected vortex beams, and their topological charges are verified.

6 Conclusion

This paper designed a bilayer metal metasurface operating in the 0.1 THz band. Based on geometric phase and far-field scattering principles, it achieved parallel generation and independent control of multi-channel vortex beams with multiple topological charges. First, a transmissive double-L-shaped metasurface unit using polyimide as the dielectric layer and copper as the metal layer was proposed. By optimizing structural parameters, high transmission efficiency near the target frequency was achieved, allowing flexible control of the cross-polarized phase response through unit rotation. Second, an interleaved phase composition method was adopted, integrating multiple independently controlled phase distributions onto a single metasurface. This effectively generated vortex beams with topological charges l=±1 and ±2, with minimal crosstalk between channels, validating the "one metasurface integrates multiple functions" design concept. Further combining phase gradient and vortex phase superposition, under LCP/RCP incidence, multi-channel vortex beams with different azimuth angles and topological charges were output. Simulation results clearly showed the helical phase and hollow intensity distribution of each channel beam. Additionally, MATLAB programming was used for automatic phase distribution calculation and CST modeling, improving design efficiency. This research provides a compact, integrated solution for efficient generation and multiplexing of vortex beams in the terahertz band, with application potential in terahertz communication, imaging, and particle manipulation. Future work could explore bandwidth expansion, efficiency improvement, and dynamic reconfigurable systems.

Reference:

1. J. W. Zhou, Y. L. Yin, J. H. Tang, Y. Xia, J. P. Yin, Frontiers of Physics 19(5), 52202 (2024).
2.  A. E. Willner, H. Huang, Y. Yan, Y. Ren, N. Ahmed, G. Xie, C. Bao, L. Li, Y. Cao, Z. Zhao, J. Wang, M. P. J. Lavery, M. Tur, S. Ramachandran, A. F. Molisch, N. Ashrafi, S. Ashrafi, Advances in Optics and Photonics 7(1), 66-106 (2015).
3. N. F. Yu, F. Capasso, Nature Materials 13(2), 139-150 (2014).
4. H. Yang, K. Ou, H. Y. Wan, Y. Q. Hu, Z. Y. Wei, H. H. Jia, X. B. Cheng, N. Liu, H. G. Duan, Materials Today 67, 424-445 (2023).
5. Y. N. Zhang, Y. G. Chen, T. Wang, Q. Zhu, M. Gu, Opto-Electronic Advances 7, 230194 (2024).
6. A. M. Berhe, K. Asham, I. Al-Ani, H. T. Hattori, A. E. Miroshnichenko, Opto-Electronic Advances 7, 230181 (2024).
7. S. X. Zhang, H. Li, C. Z. Fan, Z. C. Zeng, C. Xiong, J. Wu, Z. J. Yan, D. M. Liu, Q. Z. Sun, Opto-Electronic Advances 7, 240119 (2024).
8. Z. J. Yu, M. X. Li, Z. Y. Xing, H. Gao, Z. Y. Liu, S. L. Pu, H. Mao, H. Cai, Q. Ma, W. Q. Ren, J. Zhu, C. Zhang, Opto-Electronic Science 3, 240016 (2024).
9. R. C. Devlin, A. Ambrosio, N. A. Rubin, J. P. B. Mueller, F. Capasso,  Science 358(6365), 896-901 (2017).
10. H. Rubinsztein-Dunlop, A. Forbes, M. V. Berry, M. R. Dennis, D. L. Andrews, M. Mansuripur, C. Denz, C. Alpmann, P. Banzer, T. Bauer, E. Karimi, L. Marrucci, M. Padgett, M. Ritsch-Marte, N. M. Litchinitser, N. P. Bigelow, C. Rosales-Guzmán, A. Belmonte, J. P. Torres, T. W. Neely, M. Baker, R. Gordon, A. B. Stilgoe, J. Romero, A. G. White, R. Fickler, A. E. Willner, G. D. Xie, B. McMorran, A. M. Weiner,  Journal of Optics 19(1), 013001 (2017).
11. T. Stav, A. Faerman, E. Maguid, D. Oren, V. Kleiner, E. Hasman, M. Segev, Science 361(6407), 1101-1104 (2018).
12. N. F. Yu, P. Genevet, M. A. Kats, F. Aieta, J. P. Tetienne, F. Capasso, Z. Gaburro, Science 334(6054), 333-337 (2011).
13. Z. E. Bomzon, et al., “Pancharatnam–Berry phase in space-variant polarization-state manipulations with subwavelength gratings,” Optics Letters 26(18), 1424-1426 (2001).
14. L. L. Huang, X. Z. Chen, H. Mühlenbernd, G. X. Li, B. F. Bai, Q. F. Tan, G. F. Jin, T. Zentgraf, S. Zhang, Nano Letters 12(11), 5750-5755 (2012).
15. T. J. Cui, M. Q. Qi, X. Wan, J. Zhao, Q. Cheng, Light: Science & Applications 3(10), e218-e218 (2014).
16. Z. L. Deng, M. Jin, X. Ye, S. Wang, T. Shi, J. Deng, N. Mao, Y. Cao, B.-O. Guan, A. Alù, G. Li, and X. Li,  Advanced Functional Materials 30(21), 1910610 (2020).