1. Introduction

Vector vortex beams (VVBs), characterized by both orbital angular momentum and spatially inhomogeneous polarization, hold significant application value in cutting-edge fields such as quantum information processing, optical micromanipulation, and super-resolution imaging [1,2,3]. Their multiple degrees of freedom in phase and polarization provide a core physical foundation for high-dimensional optical communications and precise particle manipulation. However, efficient generation and flexible regulation still face technical bottlenecks: conventional spiral phase plates and spatial light modulators suffer from inherent drawbacks including low energy efficiency, limited modulation bandwidth, and insufficient system integration, making it difficult to meet the demands of complex optical field manipulation [4,5,6].

Metasurfaces, with their flexible electromagnetic responses enabled by subwavelength unit cells, have emerged as a core technological solution for vortex beam manipulation. Nevertheless, conventional optimization of unit structures relies on parameter sweeping, simulated annealing, tabu search, and other methods, which entail long design cycles, heavy computational costs, and difficulty in approaching the global optimum [7,8,9,10]. The optical field manipulation performance of metasurfaces is directly determined by the unit geometry, and the precise design of unit structures constitutes the central prerequisite for achieving high-performance vortex beam devices [11,12].

Deep learning, with its powerful nonlinear mapping and high-dimensional feature extraction capabilities, can rapidly establish the intrinsic relationship between structural parameters and optical responses, thereby providing an efficient solution pathway for inverse metasurface design [13,14,15]. Existing studies have extensively applied various deep networks to all-dielectric, metallic, and terahertz metasurface systems, demonstrating significant efficiency advantages in the design of absorbers, polarization converters, vortex beam generators, and other functional devices [16,17,18,19,20,21]. The introduction of physics-informed machine learning and cascaded network architectures has further alleviated the inherent multi-solution nature and convergence difficulties of electromagnetic inverse problems, effectively enhancing the physical reliability of inverse design results [22,23]. Meanwhile, the mature development of ReLU and LeakyReLU activation functions as well as automatic differentiation techniques has provided solid algorithmic support for the stable and efficient training of deep neural networks [24,25,26].

To address the design requirements of terahertz vector vortex beam manipulation, this paper proposes a cascaded deep learning network-based inverse design method, combined with 3D printing technology for fabricating all-dielectric alumina metasurfaces. Leveraging the synergistic regulation mechanism of propagation phase and geometric phase, this scheme enables efficient generation of vector vortex beams with different topological charges and distinct polarization characteristics, offering a feasible technological pathway for the intelligent and rapid manipulation of multi-dimensional optical fields in the terahertz regime [27,28,29,30].

2. Principles and Metasurface Design

In exploring the complex mapping between inputs and outputs, the performance of deep learning is governed by multiple factors, among which the quality of the training dataset plays a decisive role. High-quality training data endow the deep learning model with stronger learning capability, enabling it to capture effective feature representations for target prediction and thereby enhance its predictive accuracy. Moreover, a high-quality dataset helps mitigate the risk of overfitting and improves the model's generalization performance on unseen samples.

In this study, an I-shaped all-dielectric metasurface structural model is adopted. The basic constituting unit cell is shown in Fig. 1, which consists of a top I-shaped patterned layer and a substrate. The period ‘p’ is fixed at 1.6 mm, and both layers are made of alumina with a refractive index of 2.9. The substrate thickness w is 1 mm, and the height of the I-shaped pattern h is 2.5 mm. The operating frequency is chosen as 100 GHz. The top alumina pattern is composed of three rectangles, with the width of each rectangle ranging from 0.3 to 1 mm and the length exceeding 0.3 mm. For generation, Python code is employed to randomly produce coding matrices, and modeling as well as simulation are carried out via the built-in Python interface of CST according to the structural codes. The simulation frequency range is 95–105 GHz, from which the S-parameters under normal incidence in both forward and backward directions are calculated, including the reflection coefficients for both co-polarization and cross-polarization. The structures and their corresponding S-parameters are then recorded in pairs. By repeating this process iteratively, a large-scale dataset can be rapidly generated on demand.

The electromagnetic response of the metasurface is closely related to the geometric configuration of the top alumina pattern, and the arrangement matrix of the I-shaped unit cells directly determines the corresponding S-parameter characteristics. In this study, a Python random algorithm is introduced in the generation of coding matrices, enabling fully automated generation, which reduces the bias introduced by manual design while expanding the sampling space of structural parameters. Due to the extremely large number of possible permutations and combinations of the coding matrix, achieving a complete dataset would require millions of samples, which is evidently impractical under the available computational resources. Therefore, deep learning techniques are employed to compress computational costs and approach the global optimum solution.

Generally speaking, a sufficiently large data scale helps improve the generalization performance and prediction accuracy of neural networks in target scenarios; however, it also brings about a substantial increase in training time and computational consumption. Based on a comprehensive trade-off between model performance gain and computational resource investment, this study constructs a fundamental database containing 10,000 valid samples. To ensure training efficiency and the objectivity of evaluation results, this dataset is divided into three non-overlapping subsets: 70% is allocated as the training set for iterative optimization of model parameters, 20% constitutes the validation set for dynamically monitoring model performance and adjusting hyperparameters during the training phase, and the remaining 10% is reserved as the test set for independently evaluating the model's generalization performance after training. Nevertheless, the raw data must undergo normalization before being fed into the network. Specifically, the raw quantities representing specific physical parameters (S-parameters) need to be converted into data structures recognizable by the neural network, and this preprocessing step plays a critical role in ensuring the quality and reliability of the subsequent training process.

To achieve the generation of vector vortex beams, it is necessary to leverage the metasurface for flexible modulation of the phase and amplitude of the optical field. The metasurface structure employed in this chapter is composed of I-shaped unit cells. The selected material is alumina, which possesses a refractive index of 2.9 and exhibits favorable transmission efficiency as well as phase retardation at 0.1 THz. By adjusting the dimensions of each part of the unit structure, individual unit cells can be independently optimized. The unit structures obtained through training and optimization all satisfy the characteristics of half-wave plates, and their wavefront phase can cover a full range of 2π. The combination of these unit cells enables precise control over the wavefront of an incident plane wave. A synergistic approach combining geometric phase and propagation phase is adopted to regulate the optical field modulation performed by the metasurface.

The phase for generating vector vortex beams is obtained through the superposition of a series of phase components. By combining the principles of propagation phase and geometric phase, the mathematical expression can be formulated as:

(1)

(2)

(3)

where k=λ/2π is the wavenumber, λ is the wavelength, r is the polar radial coordinate, f is the focal length, Δf is the extended depth of focus, R is the radius of the metasurface, l_i and l j  are the topological charges, is the azimuthal angle, ξ is the propagation phase of the unit cell, and θ is the rotation angle of the unit cell. Equations (1) and (2) describe the phase distribution of vortex beams. Based on these equations, vortex beams with different topological charges can be generated along distinct polarization directions, or vectorially polarized vortex beams can be produced. Equation (3) provides the connection relating the vortex beam phase to the unit cell structures and their arrangement.

The overall model architecture comprises two submodules: a forward prediction network and an inverse design network. The former performs the mapping from metasurface structural parameters to optical response characteristics, while the latter carries out the inverse operation of deducing structural parameters from target optical responses. As illustrated in Fig. 2, the forward prediction network adopts a dual-input single-output architecture. Input stream 1 receives structural images that depict the spatial topological arrangement of the metasurface, whereas input stream 2 receives spectral position encodings that represent the frequency distribution information of the electromagnetic waves. Following encoding by their respective feature extraction modules, the two input streams are concatenated at the fusion layer to form a joint representation, which is subsequently passed through the response regression module to generate the predicted optical response distribution.

The core function of the forward prediction network is to establish an implicit correlation model between the metasurface topological patterns and the electromagnetic spectral responses. Upon input of a structural image, the network can output the corresponding electromagnetic response prediction within milliseconds, offering significant advantages in both timeliness and stability compared with the computational cycles of conventional full-wave simulations.

The image encoding branch comprises four groups of convolutional layers, with batch normalization and max-pooling units applied after each layer. The numbers of convolutional kernels are 32, 64, 128, and 256 successively, with a uniform kernel size of 3×3 and a stride of 1. The pooling window is set to 2×2. Through layer-wise downsampling, the spatial dimensions are compressed from 32×32 to 2×2, and subsequently flattened into a 1024-dimensional feature vector. Regularization is deployed only at the shallowest layers where information redundancy is highest; a Dropout layer with a dropout rate of 0.2 is inserted after the first convolutional group, yielding an output feature map of size 16×16×32. This random masking of a portion of channels enhances generalization capability without disrupting the core structural patterns.

The spectral embedding branch consists of three fully connected layers with neuron counts of 256, 512, and 1024, respectively, which map the discrete spectral position information into an embedding space of the same dimensionality as the image features. A Dropout layer with a dropout rate of 0.2 is inserted after the second fully connected layer in this branch, at which point the feature dimension is 512. The two branches of features are concatenated at the fusion layer to form a 2048-dimensional joint representation, which is then fed into the response regression module. This module comprises three fully connected layers with decreasing neuron counts of 1024, 512, and 1001, and ultimately outputs the electromagnetic response curve covering the entire frequency band.

All convolutional layers and fully connected layers adopt the LeakyReLU activation function with a negative slope parameter of 0.1. This activation function preserves the convergence speed advantage of ReLU while effectively alleviating the "dying neuron" problem and enhancing gradient flow during backpropagation. During the training phase, the mean squared error (MSE) is employed as the loss function, which measures model performance by computing the mean of the squared differences between predicted and true values. Its convex nature ensures stability in the optimization process and makes it well-suited for gradient descent-based algorithms to seek the global minimum. In addition, MSE is highly sensitive to small deviations, which facilitates fine fitting of spectral details, and it imposes larger penalties on outliers, thereby automatically suppressing the interference of anomalous samples during training.

The inverse design network is functionally positioned to reconstruct metasurface structural images from preset electromagnetic responses, enabling inverse deduction from performance metrics to geometric schemes and overcoming the efficiency bottleneck of conventional iterative design approaches. The network is topologically mirror-symmetric to the forward prediction network: its input layer receives a 1001-dimensional target response vector, and its output layer generates a 32×32 structural image matrix. The network architecture adopts a single-input multi-output decoder design. The front end consists of a response decoding module comprising four fully connected layers, which map the high-dimensional response vector into a feature space of the same dimensionality as the fusion layer. The back end is an image generation module consisting of four groups of deconvolution–batch normalization–upsampling units with an upsampling factor of 2, which progressively restore the spatial dimensions from 2×2 to 32×32 through layer-wise upsampling. Finally, a 1×1 convolutional layer reduces the dimension to a single channel, and a Sigmoid activation function is applied to output the binarized structural image.

The training procedure is divided into two stages. First, the forward prediction network is trained independently, and upon convergence, all its weight parameters are saved. It is then placed at the front end of the inverse design network as a fixed feature extractor, thereby constraining the solution space and alleviating the convergence difficulty arising from the many-to-one mapping from spectral responses to structural images. Both the forward prediction network and the inverse design network undergo 1,000 iterations each, and their learning dynamics are presented in Fig. 3

Fig. 3(a) illustrates the loss evolution trajectory of the forward prediction network throughout the entire training cycle. The blue curve corresponds to the MSE loss on the training set, while the red curve represents the MSE loss on the validation set. In the early stage of training, both curves decline steeply in synchronization; subsequently, the decay rate gradually slows down, and both eventually level off. At the end of training, the MSE loss of the forward prediction network on the training set drops below 0.01, and the validation set loss stabilizes at approximately 0.02. The training error and validation error maintain a reasonable gap, indicating that the network does not suffer from overfitting and exhibits good prediction stability. Fig. 3(b) presents the learning curve trend of the inverse design network. The final training loss of this network approaches 0, while the validation set loss is approximately 0.04.

To verify the prediction capability of the trained forward prediction network, we randomly selected two different metasurface structural parameter samples from the test set and compared the results obtained from CST simulation with those predicted by the model. As shown in Fig. 4, the red curves represent the spectra of the randomly selected metasurfaces from the test set obtained via CST simulation, while the blue curves denote the spectra predicted by the forward prediction network based on the metasurface structural parameters. It can be observed from the figures that the two sets of curves match well, indicating that the dataset has been adequately and effectively learned. Consequently, the forward prediction network is capable of accurately predicting the spectra of metasurfaces.

If the inverse design task were performed solely relying on a fully connected network, when the dataset contains special samples that exhibit significantly different structures yet yield similar electromagnetic responses, the network optimization process may suffer from gradient direction confusion, making stable convergence difficult to achieve. To address this many-to-one mapping dilemma, this paper adopts a cascaded architecture as a solution. The pre-trained forward prediction network is attached to the end of the inverse design network, and all its trainable parameters are frozen, allowing it to serve as a physical constraint layer.

To validate the actual design performance of the model, the target spectrum is first fed into the system, and the inverse network outputs a 32×32 structural image within milliseconds. This image is then imported into CST Microwave Studio for full-wave simulation to obtain the actual electromagnetic response, while simultaneously being fed back into the forward prediction network to obtain the theoretically predicted response. By comparing the consistency between the target spectrum (red curve) and the simulated spectrum of the inversely designed structure (blue curve), the design accuracy and physical reliability of the model are comprehensively evaluated.

Two representative sample groups are randomly selected from the test set, and the results are shown in Fig. 5. For both groups, the curves exhibit a high degree of coincidence at all frequency points, with the maximum deviation well within the acceptable range. Notably, upon verification, the structural images generated by the inverse network are found not to exist in the original dataset, confirming that the model has transcended the limitation of merely memorizing training data and possesses genuine generalization design capability. The synergistic effect of the two cascaded networks ensures a reliable solution to the inverse problem from performance metrics to geometric schemes.

3. Simulation

Based on Equations (1) and (2), the phase arrangement of the metasurface and the expected vector vortex beams are designed. In the CST software, the incident wave is set as a linearly polarized wave at a frequency of 0.1 THz, the focal length f is set to 40 mm, and Δf is set to 5 mm. The metasurface array consists of 45×45 I-shaped unit cells. In this chapter, two groups of metasurface arrays capable of generating vortex beams with distinct characteristics are designed. For the first group, under linearly polarized incidence, the left-handed circularly polarized component of the transmitted light is a vortex beam with a topological charge of +3, and the right-handed circularly polarized component is a vortex beam with a topological charge of +4. The corresponding parameters are set as l_i = 3, l_j = 4, and Δφ = 0. For the second group, also under linearly polarized incidence, the left-handed and right-handed circularly polarized components of the transmitted light are vortex beams with topological charges of +1 and −1, respectively, and the polarization state is azimuthally polarized. The corresponding parameters are set as l_i = 1, l_j = -1 and Δφ = π. Their phase arrangements are shown in Fig. 6.

First, simulation of the first metasurface array is carried out. As shown in Fig. 7, Fig. 7(a) presents the cross-sectional intensity distribution of the left-handed circularly polarized component at a distance of 40 mm from the metasurface, from which it can be observed that the intensity distribution forms a hollow ring shape, consistent with the intensity profile of a vortex beam. Fig. 7(b) shows the corresponding phase distribution at the intensity center, and through calculation, the mode purity for the topological charge l = 3 at this location is approximately 0.97, as shown in Fig. 7(c). The cross-sectional intensity distribution of the right-handed circularly polarized component at 40 mm from the metasurface is shown in Fig. 7(d), where the focused intensity also exhibits a hollow ring shape, and the phase distribution at the intensity center is presented in Fig. 7(e). The mode purity for the topological charge l = 4 is calculated to be approximately 0.98, as shown in Fig. 7(f). These simulation and calculation results are all in good agreement with the expected design targets, demonstrating that the unit structures optimized through deep learning are reliable and capable of constituting metasurface arrays for manipulating vector vortex beams at 0.1 THz.

Subsequently, to observe the polarization states of the vortex beams, the second metasurface array is also simulated, with the results shown in Fig. 8. Likewise, after linearly polarized light is incident on the metasurface, the left-handed circularly polarized component of the transmitted light is a vortex beam with a topological charge of +1, as shown in Fig. 8(a). Fig. 8(b) displays the phase distribution at this position, and the mode purity for the topological charge *l* = +1 is calculated to be approximately 0.98, as presented in Fig. 8(c). The intensity corresponding to the right-handed circularly polarized component of the transmitted light is shown in Fig. 8(d). The phase distribution and topological charge mode purity under this condition are also presented and quantitatively calculated, as shown in Fig. 8(e) and Fig. 8(f), respectively, with the mode purity for the topological charge l = −1 being approximately 0.97.

In addition, the total intensity distribution of the transmitted light is also observed, as shown in Fig. 9(a). The intensity distributions along the x-polarization and y-polarization directions are presented in Fig. 9(b) and Fig. 9(c), respectively, exhibiting a typical two-lobe structure. Fig. 9(d)–9(f) display the polarization state distributions in the vicinity of the focal plane at 40 mm, from which it can be seen that the polarization states at this distance are all azimuthally polarized.

4. Experimental Validation

In this work, 3D printing technology is selected for the fabrication of the metasurfaces. The 3D printer employed is the 3DCR-150D ceramic 3D printer manufactured by Shanghai Shuzhao Mechatronics Technology Co., Ltd., which adopts SLA stereolithography technology. The printing precision can reach 37.5 μm, with a layer thickness ranging from 0.03 to 0.1 mm and a maximum printing depth of 200 mm. The applicable printing materials include alumina, zirconia, silicon nitride ceramics, and others. The specifications of this printer fully meet the dimensional requirements of the designed unit cells.

First, the optical path for experimental testing is set up on an optical table, as shown in Fig. 10(a). The electromagnetic wave source is still the terahertz source from TeraSense. The terahertz source is connected to a horn antenna for beam emission, and the emitted wave can be regarded as a spherical wave. An aperture stop is employed to constrain the beam, followed by an off-axis parabolic mirror to correct spherical aberrations of the constrained wave. After reflection, the beam is further filtered by another aperture stop to suppress stray waves, and finally illuminates the device. This configuration ensures that the distance from the source to the sample is sufficiently large to guarantee the plane-wave characteristics of the incident wave, while also ensuring the quality of the transmitted wave with minimal stray interference. The transmitted light after conversion by the device is received by a terahertz camera, and the acquired data are processed at the terminal. Likewise, the terahertz camera is placed at the measurement distance from the sample, and the transmission performance of the entire optical path is first debugged without the sample to ensure the stability of the source output intensity. Subsequently, the performance of the sample is tested, and all measured data are normalized. Fig. 10(b) and Fig. 10(c) show the surfaces of the actual printed samples. The geometric dimensions of the unit structures are all within the design range, and the rotation angles of the individual units vary continuously, which is consistent with the design concept.

Prior to the formal testing, the system first measures the intensity distribution of the light source without the sample to ensure that the initial incident beam fully covers the imaging area. The calibration result is shown in Fig. 11(a). To control the polarization component reaching the terahertz camera, a quarter-wave plate and a linear polarizer are sequentially placed behind the sample, so that the image captured by the terahertz camera corresponds to the target optical field. First, the first group of metasurface samples is tested. By adjusting the quarter-wave plate and the linear polarizer, the left-handed and right-handed circularly polarized components of the transmitted light are obtained, with the results shown in Fig. 11(b) and Fig. 11(c), respectively. Both intensity distributions exhibit the typical hollow ring shape characteristic of vortex beams. Subsequently, the second group of metasurface samples is tested. For this group, the terahertz camera is placed directly behind the sample for data acquisition, yielding the total transmitted intensity distribution, as shown in Fig. 11(d). A linear polarizer is then placed behind the sample, and by adjusting its orientation, the intensity distributions along the x-polarization and y-polarization directions are obtained, as shown in Fig. 11(e) and Fig. 11(f), respectively. These intensity distributions display a two-lobe structure, which is consistent with the simulation results described above.

The experimental test results show that, under linearly polarized incidence, the intensity distributions of the transmitted optical fields from the two groups of metasurface samples are in excellent agreement with the simulation results. This outcome verifies the reliability of the unit structures designed through deep learning in this chapter, demonstrating that the metasurfaces composed of these units can achieve the intended design objectives. The minor discrepancies between the experimental measurements and simulation data are mainly attributable to the fact that the actual light source in the experiment is not an ideal Gaussian beam, as well as to environmental noise interference and systematic errors introduced during the calibration of optical components.

5. Conclusions

In this chapter, deep learning is employed to optimize and identify suitable unit structures, and 3D printing technology is used to fabricate metasurfaces for generating vector vortex beams. Alumina, which exhibits favorable refractive index and transmission efficiency in the terahertz band, is selected as the material, and two groups of transmission-type metasurfaces are designed and fabricated. The phase distribution of the metasurface units is designed based on the combined principles of propagation phase and geometric phase. The intensity distribution, phase distribution, and mode purity corresponding to the topological charges of the transmitted optical field are calculated through simulations, successfully achieving vortex beams carrying different topological charges as well as azimuthally polarized states. The feasibility of the design scheme is validated through experimental measurements, and the experimental results show that the intensity distributions are consistent with the simulation findings. This fully demonstrates the feasibility of utilizing deep learning for unit structure optimization and provides important technical insights and references for the future realization of flexible and versatile light field manipulation.

Reference:

1. Y. Liu, X. Zhang, Applied Physics Letters 103, 141101 (2013).
2. G. Wang, B. Hu, M. I. Khan, et al., Frontiers in Physics 9, 751026 (2021).
3. K. Qu, K. Chen, Q. Hu, et al., Advanced Photonics Nexus 2, 016010 (2023).
4. F. Yang, C. Zheng, J. Li, et al., Optics & Laser Technology 165, 109563 (2023).
5. X. Chen, W. Su, Z. Geng, et al., Optik 267, 169701 (2022).
6. Z. Li, C. S. Prasad, X. Wang, et al., Nanophotonics 12, 3721–3727 (2023).
7. F. Glover, ORSA Journal on Computing 1, 190–206 (1989).
8. P. R. Wiecha, O. L. Muskens, Nano Letters 20, 329–338 (2019).
9. Z. Sun, B. Xu, F. Jin, et al., IEEE Journal of Selected Topics in Quantum Electronics 28, 1–9 (2021).
10. S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, Science 220, 671–680 (1983).
11. N. F. Yu, P. Genevet, M. A. Kats, et al., Science 334, 333–337 (2011).
12. D. R. Smith, D. C. Vier, W. J. Padilla, et al., Physical Review E 71, 036609 (2005).
13. Y. Bengio, Y. LeCun, G. Hinton, Communications of the ACM 64, 58–65 (2021).
14. A. Krizhevsky, I. Sutskever, G. E. Hinton, Advances in Neural Information Processing Systems 25, 1097–1105 (2012).
15. I. Sajedian, J. Kim, J. Rho, Microsystems & Nanoengineering 5, 27 (2019).
16. S. An, C. Fowler, B. Zheng, et al., ACS Photonics 6, 3196–3207 (2019).
17. Z. Liu, D. Zhu, S. P. Rodrigues, et al., Nano Letters 18, 6570–6576 (2018).
18. W. Ma, F. Cheng, Y. Xu, et al., Advanced Materials 31, 1901111 (2019).
19. Y. Lv, S. Liu, J. Tian, et al., Results in Physics 53, 107027 (2023).
20. W. Xi, Y. J. Lee, S. Yu, et al., Nature Communications 14, 4694 (2023).
21. B. Zheng, H. Zhang, in 2020 IEEE International Symposium on Antennas and Propagation and North American Radio Science Meeting (IEEE, 2020), pp. 2049–2050.
22. R. Pestourie, Y. Mroueh, T. V. Nguyen, et al., npj Computational Materials 6, 164 (2020).
23. L. Wang, J. Dong, W. Zhang, et al., Nanomaterials 13, 1030 (2023).
24. G. E. Karniadakis, I. G. Kevrekidis, L. Lu, et al., Nature Reviews Physics 3, 422–440 (2021).
25. D. P. Kingma, M. Welling, arXiv:1312.6114 (2013).
26. V. Nair, G. E. Hinton, in Proceedings of the 27th International Conference on Machine Learning (Omnipress, 2010), pp. 807–814.
27. A. G. Baydin, B. A. Pearlmutter, A. A. Radul, et al., Journal of Machine Learning Research 18, 5595–5637 (2017).
28. A. L. Maas, A. Y. Hannun, A. Y. Ng, in Proceedings of the 30th International Conference on Machine Learning (ICML, 2013), p. 3.
29. Y. Teng, C. Li, S. Li, et al., Optics & Laser Technology 160, 109058 (2023).
30. D. Gostimirovic, W. N. Ye, IEEE Journal of Selected Topics in Quantum Electronics 25, 1–5 (2019).