In recent years, vector vortex beams (VVBs) have attracted widespread attention in fields such as quantum information [1,2], optical trapping [3,4], and particle manipulation [5,6] due to their unique spatial intensity distribution and phase structure. VVBs not only carry orbital angular momentum (OAM) but also exhibit spatially varying polarization states [7–11], a characteristic that makes them highly valuable in scientific research and technological applications [12–14]. The introduction of OAM enables precise manipulation of microscopic particles, offering new possibilities for the development of optical tweezers [15]. Additionally, the multi-degree-of-freedom nature of VVBs holds potential applications in quantum information processing, laying the foundation for more efficient quantum communication and quantum computing [16–20]. In the fields of optical imaging and microscopy, the unique spatial structure of VVBs can significantly enhance microscope resolution, enabling super-resolution imaging [21]. Despite significant progress in VVB research, challenges remain in their generation and manipulation. Traditional methods such as spiral phase plates and spatial light modulators (SLMs) face limitations in efficiency and flexibility [22–26], necessitating further improvements.

Metasurfaces, composed of subwavelength structural units, are two-dimensional planar structures [27–30] that have emerged as a novel class of optical materials. They exhibit immense potential in manipulating wavefronts, polarization states, and phases [31–40]. These units can be designed to achieve precise control over light, thereby generating beams with specific properties [41–43]. In particular, all-dielectric metasurfaces, utilizing materials like aluminum oxide, feature low absorption losses and high refractive indices, making them suitable for high-performance optical devices [44]. The advancement of 3D printing technology has further propelled the fabrication of metasurfaces, enabling high-precision manufacturing of complex nanostructures.

This study proposes an all-dielectric aluminum oxide metasurface fabricated using 3D printing technology for efficiently generating longitudinally varying vortex beams. Aluminum oxide was chosen due to its relatively low cost, excellent optical properties, and high refractive index, making it an ideal material for optical components. Additionally, 3D printing technology, being mature and cost-effective, enables the precise fabrication of complex structures, providing significant convenience for designing and manufacturing metasurfaces with specific functionalities. We utilized stereo lithography appearance (SLA) technology to precisely control the micro-nano structures of the metasurface, thereby achieving customized modulation of the phase of the transmitted beam and generating vortex beams with desired characteristics. Compared to traditional 3D printing techniques, SLA offers advantages such as high precision, fast speed, and low cost.

Based on the principles of Pancharatnam–Berry (PB) phase [45,46] and propagation phase [47,48], we designed an all-dielectric aluminum oxide metasurface operating at a frequency of 0.1 THz to efficiently generate longitudinally varying vector vortex beams with high precision. By leveraging the two degrees of freedom provided by the propagation phase and PB phase, we achieved flexible control over the field distribution of orthogonal circular polarization states by superimposing two vector vortex beams with different focal lengths in different polarization directions, as illustrated in Fig. 1(a). The PB phase refers to the phenomenon where, when illuminated by circularly polarized light, the rotation angle (counterclockwise) of the unit structure, denoted as α, induces a polarization change in the anisotropic medium, causing the cross-polarized component to carry an additional phase factor of exp(2iα). This allows phase modulation by rotating the unit structures. On the other hand, the propagation phase is achieved by designing the geometric shape and size of the unit structures, enabling the accumulation of phase differences as light waves pass through them, thus achieving fine control over the phase of the light waves.

We utilized 3D-printable aluminum oxide (Al2O3) to construct the metasurface. Aluminum oxide has a high refractive index of 2.9 at 0.1 THz, and its low-loss characteristics, excellent stability, and high transmission in the terahertz band make it an ideal material for phase modulation. These properties enable aluminum oxide to effectively achieve phase modulation in the terahertz frequency range. The unit structure we employed is a rectangular pillar with a periodicity of Px=Py=1.6mm and a height of h=2.5mm, as shown in Fig. 1(b). The substrate has a thickness of w=1mm. The unit structure was optimized by independently adjusting the dimensions of the long and short axes (x and y).

Considering that the metasurface functions similarly to a phase retarder, when the long and short axes of the designed structure align with the Cartesian coordinate system, the transmission matrix can be expressed as T=[Ayeiφy00Axeiφx],(1)where Ax and φx represent the amplitude and phase of the incident light polarized along the x-axis, while Ay and φy denote the amplitude and phase of the incident light polarized along the y-axis. When the unit structure is rotated by an angle θ, the long and short axes no longer align with the x- and y-axes. In this case, the transmission matrix transforms into T=[cosθ−sinθsinθcosθ][A2eiφ200A1eiφ1][cosθsinθ−sinθcosθ].(2)

Since the long and short axes no longer align with the x- and y-axes, A1 and φ1 represent the amplitude and phase of the incident light polarized along the x-axis when the structure is not rotated, while A2 and φ2 represent the amplitude and phase of the incident light polarized along the y-axis when the structure is not rotated. Considering a special case where A1≈A2≈A and |φ1−φ2|≈π, according to Euler’s formula, we have exp(iφ1)=−exp(iφ2). Therefore, Eq. (2) can be simplified as $$\begin{gathered} T=Aeiφ2[cosθ−sinθsinθcosθ][100−1][cosθsinθ−sinθcosθ] \\ =Aeiφ2[cos2θsin2θsin2θ−cos2θ]. \end{gathered}$$(3)

From Eq. (3), it can be concluded that under the conditions A1≈A2≈A and |φ1−φ2|≈π, modifying the dimensions of the long or short axis of the unit structure or the rotation angle can alter the transmission matrix, thereby affecting the outgoing light. Therefore, when searching for unit structures, it is essential to ensure high transmission efficiency and satisfy |φ1−φ2|≈π, effectively functioning as a half-wave plate. This satisfies the Pancharatnam–Berry (PB) phase theory, where the cross-polarized component of the transmitted beam carries an additional phase factor of exp(2iθ) when the unit structure is rotated by θ. Furthermore, φ1 and φ2 must cover a range of 2π to enable wavefront manipulation, adhering to the transmission phase theory.

We explored the dimensions of the unit structure within the range of 0.1–1.4 mm for the long and short axes. Using CST Studio Suite, we simulated the unit structures at 0.1 THz and identified 16 unit structures that meet the requirements; detailed parameters can be found in Appendix A. Their transmission efficiencies Txx and Tyy are both above 0.85, the phase differences |φx−φy| are approximately 180°, and they cover a full 360° phase range, as shown in Fig. 1(d). To verify that these units satisfy the PB phase theory, we adjusted the incident wave to circularly polarized light and rotated the unit structures by α=45°, 90°, 135°, and 180°, as illustrated in Fig. 1(c). The simulation results, shown in Figs. 1(e) and 1(f) (the “serial number” represents the phase of the unit structure arranged from small to large), demonstrate that the phase changes are approximately 2α, and a rotation of 180° covers a full 360° phase shift. Thus, the dimensions of these unit structures are reliable.

Among the various vortex beams, the Laguerre–Gaussian (LG) beam is a typical type of vortex beam. Due to its well-defined mathematical expression and physical properties, it is widely used in various optical experiments and theoretical research. The Laguerre–Gaussian beam has a spiral phase structure in space and can carry a well-defined orbital angular momentum (OAM), with a phase distribution of exp(imφ). In cylindrical coordinates, the electric field expression of an LG beam propagating along the axial direction can be expressed as [49] Epm(ρ,φ,ζ)=2p!π(p+|m|)!1ω(ζ)[2ρω(ζ)]|l|exp[−ρ2ω2(ζ)]⁢Mp|m|[2ρ2ω2(ζ)]exp(imφ)exp[ikρ2ζ2(ζ2+ζR2)−i(2p+|m|+1)arctan(ζζR)],(4)ω(ζ)=ζ01+(ζζR)2,(5)ζR=πω02λ,(6)where m represents the topological charge, which is the numerical value describing the OAM quantum state; p denotes the radial mode number, used to define the radial mode of the beam; while ρ and φ correspond to the radial and angular coordinates in cylindrical coordinates, respectively. Additionally, ζ describes the distance along the propagation direction of the beam, ω(ζ) is the beam radius, Mp|m| is the generalized Laguerre polynomial related to p and m, the Rayleigh range is denoted by ζR, and the wave number k is defined as 2π/λ. For circularly polarized LG beams, it can be expressed as Eτ=±1L/R=Epm·[1±i]=Bpmexp(imφ)·[1±i],(7)where τ=±1 indicates the states of left-handed circularly polarized (LCP) light and right-handed circularly polarized (RCP) light, and Bpm is the coefficient of the LG beam in Eq. (4). When the radial mode number p=0 and two LG beams with opposite topological charges are superimposed, according to Eq. (7), we obtain Esum=Bmexp(im1φ)·[1i]+B−mexp(−imφ)·[1−i]=2Bm[cos(mφ)−sin(mφ)]=ERad.(8)

At this point, such superposition of two LG beams constitutes a radially polarized vortex beam. Furthermore, if the phase difference between these two LG beams is π, then Eq. (8) can be rewritten as Esum=Bmexp(im1φ)·exp(iπ)·[1i]+B−mexp(−im2φ)·[1−i]=−2iBm[sin(mφ)cos(mφ)]=EAz.(9)

In this case, the superposition of the two LG beams results in an azimuthally polarized vector vortex beam. To generate longitudinally varying vortex beams, the phase distribution of the circularly polarized channels needs to be defined as [35] $$\begin{gathered} ψjL/R(ρ)=k·ρ22(fj+Δf·ρ2/R2)+mjL/Rφ, \\ ΦL/R(ρ)=arg{∑j=1Nexp[iψjL/R(ρ)]}, \end{gathered}$$(10)where Δf is the extended depth of focus, fj is the focal length, and R is the radius of the metasurface. If generating longitudinally varying vector vortex beams, Eq. (10) can be rewritten as $$\begin{gathered} ψjL(ρ)=k·ρ22(fj+Δf·ρ2/R2)+mjφ+Δφ, \\ ψjR(ρ)=k·ρ22(fj+Δf·ρ2/R2)−mjφ, \\ ΦL/R(ρ)=arg{∑j=1Nexp[iψjL/R(ρ)]}, \end{gathered}$$(11)where Δφ is the phase difference between the two circularly polarized channels. Furthermore, as derived above, when Δφ=0, the generated vector vortex beam is radially polarized, and when Δφ=π, the generated vector vortex beam is azimuthally polarized.

A linearly polarized plane wave can be decomposed into two circularly polarized components with opposite handedness, namely, LCP and RCP, considering its incidence: $$\begin{gathered} EinL=22ELeiδL[1i], \\ EinR=22EReiδR[1−i]. \end{gathered}$$(12)

According to Eqs. (3) and (12), the outgoing light after passing through the metasurface can be expressed as $$\begin{gathered} EoutLR=T·EinL=22AELeiδLeiφ2ei2θ[1−i]=22AELeiδLei(φ2+2θ)[1−i], \\ EoutRL=T·EinR=22AEReiδReiφ2e−i2θ[1i]=22AEReiδRei(φ2−2θ)[1i], \end{gathered}$$(13)where EoutLR represents the right-handed circularly polarized component transmitted after LCP incidence, EoutRL represents the left-handed circularly polarized component transmitted after RCP incidence, φ2 is the phase shift of the unit structure under y-axis polarized incidence, and ±2θ is the geometric phase.

First, we designed a metasurface capable of realizing longitudinally varying vortex beams. The metasurface consists of an array of 45×45 unit structures, operating at a frequency of 0.1 THz. Its function is to generate vortex beams with topological charges longitudinally varying from +3 to −3 for LCP incidence and from +4 to −4 for RCP incidence. The phase distribution of the metasurface can be given by Eq. (10). For the left-handed circularly polarized component ΦL(ρ) and the right-handed circularly polarized component ΦR(ρ), two vortex beams with different topological charges can be superimposed at different distances, namely, ψ1L(ρ), ψ2L(ρ), ψ1R(ρ), and ψ2R(ρ), where k=2π/λ. The focal lengths can be set as f1=22mm, f2=42mm, and Δf=5mm. The final phase distribution of the metasurface is shown in Fig. 2. Combining Eq. (13), the phase shift and rotation angle of the unit structure are given by Eq. (14): $$\begin{gathered} φ2−2θ=ΦL(ρ), \\ φ2+2θ=ΦR(ρ). \end{gathered}$$(14)

Based on the phase distribution shown in Fig. 2, we designed the 45×45 units metasurface. In CST Studio Suite software, the time-domain solver was employed for simulation calculations. The incident waves were LCP and RCP at 0.1 THz. The simulation results provided the electric field distribution on the xoy plane at different distances. Figures 3(a) and 3(b) respectively illustrate the intensity distribution and phase distribution of the vortex beam under LCP incidence as the distance varies. Between 22 mm and 42 mm, the topological charge evolves from m=+3 to m=−3. This is due to the focusing of vortex beams with different topological charges at different distances, and the transition is gradual with intermediate states. Similarly, Figs. 3(c) and 3(d) show the intensity distribution and phase distribution of the vortex beam under RCP incidence. Between 22 mm and 42 mm, the topological charge evolves from m=+4 to m=−4.

Additionally, we can employ Fourier transform analysis to quantitatively evaluate the mode purity of the topological charge in the vortex beams: $$\begin{gathered} E(φ)=∑mAmeimφ, \\ Am=12π∫02πE(φ)exp(−imφ)dφ, \\ Pm=|Am|2∑m=−∞m=∞|Am|2, \end{gathered}$$(15)where E(φ) is the angular distribution function of the optical field, Am represents the Fourier coefficients of E(φ), and Pm denotes the energy weight of the OAM mode in the total energy, i.e., the mode purity.

Figure 3(e) illustrates the variation of mode purity with distance for the vortex beams generated under LCP and RCP incidence. Clearly, at propagation distances of 22 mm and 42 mm, for LCP incidence, the mode purity is primarily composed of +3 and −3. As the propagation distance increases, the mode purity gradually shifts from being dominated by +3 to being dominated by −3. Similarly, for RCP incidence, the mode purity is mainly composed of +4 and −4, and as the propagation distance increases, it transitions from being dominated by +4 to being dominated by −4. This quantitative analysis confirms that the simulation results align with the expected outcomes.

Furthermore, by combining Eqs. (8) and (9), a pair of orthogonal circularly polarized light beams can be designed to generate vortex beams with opposite topological charges and a phase difference at the near focal plane, and vortex beams with opposite topological charges at the far focal plane. This enables the realization of longitudinally varying vector vortex beams, which exhibit azimuthal polarization at the near focal plane and radial polarization at the far focal plane. The phase distribution of the metasurface can be expressed by Eq. (11), with the focal lengths f1, f2, and Δf kept constant, mj=1, and Δφ=π. The phase shifts and rotation angles of the specific unit structures are still determined by Eq. (14).

Figure 4 demonstrates the transformation of the polarization state from azimuthal to radial polarization as the propagation distance increases from 22 to 42 mm under y-linearly polarized light incidence. Figures 4(a) and 4(b) illustrate the evolution of the electric field intensity distributions in the x- and y-polarization directions at different distances. In the x-polarization mode, the initial double-lobe structure at 22 mm is symmetrically distributed along the y-axis. As the propagation distance increases, it undergoes continuous rotation, completing a 90° spatial orientation transformation at 42 mm to form a double-lobe structure symmetric along the x-axis. Conversely, the y-polarization channel exhibits an orthogonal symmetric evolution, with its double-lobe structure transitioning from an initial x-axis symmetry to a y-axis symmetry. Figure 4(c) shows the specific distribution of the polarization states, where the azimuthal polarization at 22 mm evolves into radial polarization at 42 mm, with transitional characteristics observed near both 22 mm and 42 mm. This fully validates the feasibility of longitudinally varying vector vortex beams.

To further validate that the proposed metasurface can achieve the evolution of the polarization state of vector vortex beams, we designed two sets of metasurface configurations. These configurations are intended to realize the evolution from azimuthal polarization to second-order radial polarization and from radial polarization to second-order radial polarization, respectively.

For the evolution from azimuthal polarization to second-order radial polarization, as indicated by Eqs. (8) and (9), it is necessary to configure two circularly polarized beams with opposite topological charges and a phase difference π at the near focal plane, and beams with opposite topological charges and m=±2 at the far focal plane, as described in Eq. (16): $$\begin{gathered} ψ1L(ρ)=k·ρ22(fj+Δf·ρ2/R2)+φ+π,ψ2L(ρ)=k·ρ22(fj+Δf·ρ2/R2)+2φ, \\ ψ1R(ρ)=k·ρ22(fj+Δf·ρ2/R2)−φ,ψ2R(ρ)=k·ρ22(fj+Δf·ρ2/R2)−2φ. \end{gathered}$$(16)

Figures 5(a), 5(b), and 5(e) demonstrate the intensity distributions and polarization state distributions in different polarization directions as the output light transitions from azimuthal polarization to second-order radial polarization. The intensity distribution of the transmitted x-polarization evolves from an initial double-lobe structure symmetric along the y-axis at 22 mm to a composite petal pattern with complementary x/y-axis distributions at 42 mm. Correspondingly, the y-polarization channel exhibits an orthogonal symmetric evolution, transitioning from an initial double-lobe structure symmetric along the x-axis to a composite petal pattern with complementary x/y-axis distributions. This complementary evolution characteristic of the orthogonal polarization channels is directly validated in the polarization state distribution shown in Fig. 5(e), where the initial axisymmetric feature of azimuthal polarization evolves into a second-order radial polarization with four-quadrant polarization. The polarization vectors in adjacent quadrants exhibit a π/2 phase difference, perfectly corresponding to the complementary petal characteristics of the intensity distribution.

Similarly, Figs. 5(c), (d), and 5(f) illustrate the intensity distributions and polarization state distributions in different polarization directions as the output light transitions from radial polarization to second-order radial polarization. The intensity distribution of the transmitted x-polarization evolves from an initial double-lobe structure symmetric along the x-axis to a composite petal pattern with complementary x/y-axis distributions. Meanwhile, the y-polarization channel continues to exhibit an orthogonal symmetric evolution. This polarization order doubling effect is perfectly represented in the polarization distribution shown in Fig. 5(f). Both sets of simulation results precisely align with the theoretical design of the metasurface, fully validating the effectiveness of the longitudinally varying vector vortex beam design.

We have chosen 3D printing technology to fabricate the metasurface, using the 3DCR-150D ceramic 3D printer from Shanghai Digital Manufacturing Technology Co., Ltd. This printer employs SLA technology, achieving a printing accuracy of up to 37.5 μm, with a layer thickness ranging from 0.03 to 0.1 mm and a maximum printing depth of 200 mm. It is compatible with materials such as alumina, zirconia, and silicon nitride ceramics. The printer’s specifications fully meet the requirements of our designed metasurface unit dimensions.

We constructed the experimental test optical setup on an optical platform, as shown in Fig. 6(a). The electromagnetic wave source used was a terahertz source from TeraSense Company, with a frequency range of 0.1–0.3 THz. The terahertz source was connected to a horn antenna (model Anteral SGH-26-WR10) for beam emission, with the source operating in a point-source emission mode, emitting spherical waves. An aperture diaphragm was used to constrain the beam, followed by an off-axis parabolic mirror (model MPD369-M01 from Thorlabs) to correct spherical aberrations in the constrained wave. The reflected beam was further filtered by another aperture diaphragm to remove stray waves before illuminating the device. This setup ensured a sufficiently large distance between the source and the sample, maintaining the plane-wave characteristics of the incident wave while minimizing stray wave interference. The transformed output light from the device was received by a terahertz camera (model Tera-1024 from TeraSense), and the received data was processed at the terminal. The terahertz camera operates in a frequency range of 50 GHz to 0.7 THz, with a maximum high-speed image acquisition rate of 50 frames per second. The camera was placed at the measurement sample distance, and the transmission performance of the entire optical path was first calibrated without the sample to ensure the stability of the source output intensity. Subsequently, the sample’s performance was tested, and all test results were normalized.

Figure 6(a₁) illustrates the distribution of a portion of the metasurface array units, where the unit structures exhibit non-uniform spatial characteristics with varying geometric dimensions and rotation angles, aligning with the design concept. Figure 6(a₂) demonstrates the excellent electromagnetic properties of the 3D-printed alumina material at a working frequency of 0.1 THz. The measured values of its refractive index (n=2.9) and extinction coefficient (k=0.02) closely match the simulation model. This high refractive index and low-loss material characteristic provide a physical foundation for the efficient modulation of the metasurface device.

Based on the aforementioned design parameters, we fabricated two sets of metasurface devices using 3D printing technology: the first set (Device I) achieves the transformation from azimuthal polarization to radial polarization, while the second set (Device II) realizes the evolution from azimuthal polarization to second-order radial polarization. The surface morphology of the devices was measured by Panorama Optical Instruments (Hangzhou) Co., Ltd., using the Olympus DSX2000 3D digital microscope from Japan. Figures 6(b) and 6(c) systematically present the surface morphology characteristics of the two devices: Figs. 6(b₁) and 6(c₁) display the surface morphology images obtained through laser scanning 3D imaging technology, while Figs. 6(b₂) and 6(c₂) provide top views, with red arrows indicating the measurement direction of surface height. Insets show schematic diagrams of the metamaterial structure height. The specific height measurement results are shown in Figs. 6(b₃) and 6(c₃), with an average height of 2.35 mm (standard deviation σ=0.05mm) for the metamaterial, deviating within 6% from the design value of 2.5 mm, confirming the precision of the fabrication process and the feasibility of the metasurface devices.

First, we recorded the intensity distribution of the light source without the metasurface device. By adjusting the experimental components, the initial light spot was made to precisely cover the imaging area, as shown on the left side of Fig. 7(a). To accurately characterize the vector modulation properties of the device, a linear polarizer (PW010-025-075) was added in front of the terahertz camera to separate different polarization components. For the test results of Device I, as shown in Fig. 7(b), under y-linearly polarized light incidence, the transmitted x-polarization channel exhibited a typical y-axis symmetric double-lobe structure at 22 mm, which gradually evolved into an x-axis symmetric double-lobe structure at 42 mm. The corresponding transmitted y-polarization channel displayed an orthogonal symmetric evolution pattern, as shown in Fig. 7(c). For the test results of Device II, as shown in Fig. 7(d), under the same incident conditions, the transmitted x-polarization channel evolved from an initial y-axis symmetric double-lobe structure into a four-lobe mode with complementary x/y characteristics, while the y-polarization channel simultaneously exhibited an orthogonal symmetric evolution pattern, as shown in Fig. 7(e). As a supplement, Fig. 7(a) also illustrates the vector vortex characteristics of Device I and Device II at 22 mm and 42 mm; specific polarization calculations are provided in Appendix B. For Device I, the vortex characteristic at 22 mm was azimuthal polarization, while at 42 mm, it was radial polarization. For Device II, the vortex characteristic at 22 mm was azimuthal polarization, while at 42 mm, it was second-order radial polarization. These experimental results are largely consistent with the aforementioned simulation results, sufficiently proving that our designed structure achieves the intended functionality. The overall loss of the metasurface and quantitative comparisons are detailed in Appendix C and Appendix D.

In summary, in this study, we successfully achieved the generation of longitudinally varying vector vortex beams on an all-dielectric metasurface platform using 3D printing technology. By effectively combining the dual degrees of freedom of transmission phase and geometric phase, we integrated two vortex beams with different focal lengths and long focal depth characteristics into a pair of orthogonal circular polarization states. As a proof of principle, we designed and simulated a series of metasurface models, demonstrating the evolution of topological charges and vector polarization states of the vector vortex beams.

The experimental characterization data quantitatively matched the simulation model (error < 6%), which not only validated our theoretical expectations but also highlighted the unique advantages of 3D printing technology in the fabrication of complex micro-nano structures. However, despite the significant advantages of 3D printing technology in metasurface fabrication, it still faces certain challenges and limitations. For example, the resolution limitations of 3D printing may affect the realization of fine structures in metasurfaces, while material uniformity and printing precision could also impact device performance. Despite these challenges, 3D-printed metasurfaces hold promise for broader applications in fields such as polarization optics and optical field modulation. Our research also provides new insights into phase modulation in the terahertz frequency range.