Skip to content
Frontiers
Search

Original Research

Frontiers in Physics ·

Volume 14 - 2026 | https://doi.org/10.3389/fphy.2026.1916847

Inverse design of multi-mode OAM-generating pixelated THz metasurfaces

  • Ruilu Huang1
  • Jiusheng Li1*corresponding author
  • Guo-Hua Qiu1
  1. Centre for THz Research, China Jiliang University, Hangzhou, China
Article metrics0 viewsView details
PDFCiteMetrics

Abstract

Traditional metasurface design relies on physical priors and parameter scanning, which is time-consuming and prone to getting stuck in local optima. In this article, we proposed a pixelated metasurface inverse design method, which discretizes the metasurface elements into a 16 × 16 binary pixel array, with each pixel independently adjustable, expanding the design space. We constructed a deep neural network framework consisting of a forward prediction residual convolutional network and a tandem generative inverse network: The forward network takes the element structure image as input and quickly predicts its complex reflection coefficient spectrum. The average mean square error (MSE) of the test set is 0.0017. The inverse network is implemented by concatenating pre-trained forward networks without the need for real geometric labels, with an average MSE of 0.0131 on the test set. Using the trained model, we designed a metasurface that generates orbital angular momentum (OAM) with topological charges of ±1 and ±2 in the 1.5-1.6 THz frequency range. The full wave simulation results confirm that the designed metasurface can achieve spatial separation and independent control of multi-topology OAM beams, and the far-field amplitude and phase distribution are consistent with theoretical expectations. This method significantly improves the design efficiency and control accuracy of multifunctional metasurfaces, providing a new approach for on-demand customization of terahertz wave devices.

Keywords

1 Introduction

Metasurfaces are artificial two-dimensional electromagnetic materials composed of sub wavelength unit arrays, which can flexibly control the amplitude, phase, and polarization state of electromagnetic waves. They have shown broad application prospects in fields such as sensors, absorbers, and vortex beam generation [1–6]. In recent years, due to the ability of deep learning algorithms to quickly establish nonlinear mappings between structural parameters and electromagnetic responses, high-precision prediction results can be obtained at a lower computational cost. They have been introduced as inverse design tools in metasurface research, effectively breaking through traditional design bottlenecks [7–9]. For example, Niu et al. [10] used a tandem neural network to achieve inverse design of 1-bit encoding elements and metasurface array arrangements with 10 dB in ultra-wideband radar cross section (RCS) in 3.7–18.7 GHz. Qi et al. [11] proposed a inverse topology design method based on multi-layer perceptrons, which achieves wideband absorption with reflection attenuation ≥10 dB in 1.33–7.31 GHz band. Under incident angles less than 45°, the structure still maintains wideband absorption characteristics under TE and TM waves incidence. Ding et al. [12] used deep learning methods to carry out inverse design of a double-layer prism metasurface absorber, achieving an absorption rate greater than 97.75% in the wavelength range of 280–4000 nm. Moreover, when the incident angle is greater than 70°, the absorption rate is still greater than 92%, and the parameter design error is only 0.06%. Yan et al. [13] presented a deep learning driven polarization multiplexed metasurface intelligent inverse design framework to achieve polarization modulation switching between O-band and C-band transmission peaks. More recently, the deep learning inverse design method for metasurfaces has been extended to the terahertz frequency band. Wu et al. [14] employed a CNN Transformer hybrid network and particle swarm optimization algorithm to inverse design a metasurface structure, obtaining high-precision predictions with phase errors less than 1.03° and amplitude errors less than 2.8 × 10−3. Within a bandwidth of 0.25–0.55 THz, the phase modulation range reached 310° and the reflection amplitude was greater than 0.9, achieving directional beamforming. Chen et al. [15] constructed a inverse design model based on conditional generative networks to design chiral metasurfaces. OAM beam focusing was achieved through PB phase modulation at 1.126 THz, and the circular dichroism (CD) reached 0.471, significantly improving the design speed and accuracy of terahertz metasurfaces. Although the above-mentioned deep learning based auxiliary design methods have demonstrated excellent performance in various metasurface applications, they are all based on a common premise: metasurface elements are described by a small number of continuous geometric parameters extracted from predefined basic shapes such as rectangles, ellipses, and crosses. In such methods, the design freedom of metasurface structures is inherently limited within the low dimensional geometric space generated by the selected shape family, and cannot escape from the preset template. The overall structural design freedom is low, and the design space is limited [16–18]. To address this issue, we proposed a pixelated metasurface inverse design method for transmission phase modulation, which discretizes the metasurface elements into a 16 × 16 binary pixel array, with each pixel independently defined as metal or air. In order to construct usable datasets with limited computing resources, the D2 symmetry group constraint is introduced to expand the design space of metasurfaces to 249 possible structures. On this basis, a deep learning framework is proposed. By concatenating pre-trained forward prediction networks with frozen parameters, the inverse generation network can infer the unit geometry of metasurfaces without any real structural labels, relying solely on electromagnetic response targets. Using this framework, the successfully designed metasurface can simultaneously generate and spatially separate the orbital angular momentum with topological charges of −2, −1, +1, +2 in the frequency range of 1.5–1.6 THz. The full wave simulation results verify that the metasurface has excellent wavefront control ability. Compared with traditional metasurface design methods, this method significantly accelerates the design speed of metasurface structure, reduces trial and error time, breaks the shape constraints of traditional metasurfaces on electromagnetic waves, and opens up a new path towards phase controlled metasurface design.

2 Designing metasurfaces using deep learning methods

The designed metasurface consists of a metal pattern layer, a polyimide dielectric layer, and a metal plate bottom layer from top to bottom, as shown in Figure 1. The top pattern layer can be represented by a 16 × 16 pixels binary image, where light blue pixels represent gold (defined as “1″) and white pixels represent air (defined as “0″). Any geometric structural element design can be represented by a 16 × 16 matrix with different combinations of 0 and 1. To prevent the gold patterns in adjacent units from connecting to each other, fix the pixels at the boundaries as unfilled. The period of the metasurface element is P = 48 μm, the length of each pixel block is P1 = 3 μm, the metal thickness is set to 0.1 μm, and the dielectric layer thickness is 30 μm. Due to the fixed boundary pixels, the actual designed effective pixel area is 14 × 14, theoretically allowing for 2196 possible metasurface element patterns. Due to the limitations of current computing power, the D2 symmetry group method has been introduced in metasurface design to reduce design complexity. Firstly, the pixel pattern of the first quadrant is generated, and then it sequentially about the y-axis and x-axis is mirrored to generate complete metasurface elements. Using this method, randomly 12,000 sets of metasurface element patterns are generated for simulation.

The constructed pixelization metasurface element.
Figure 1. The constructed pixelization metasurface element.

Figure 2 shows the process of designing a pixelated metasurface based on deep learning assistance. Firstly, we use Python to generate a dataset containing 12,000 sets of metasurface unit structure images. Then, through the joint simulation of PyCharm and CST Studio Suite 2024, S-parameter data corresponding to the unit structure image is collected, and the real and imaginary parts of the S11 parameters are extracted as the output of the forward prediction neural network and the input dataset of the generative inverse network. In this article, a frequency domain solver is used, with Unit Cell boundaries in the x/y direction and Open (add space) boundaries in the z direction. We set the spectral range of 1.0–3.0 THz and discretize it into 201 frequency points as input or output data for the network. And then, we divide the data into training, validation, and testing sets in a 7:2:1 ratio, and build forward predictive neural network and generative inverse network based on PyTorch framework for designing metasurfaces and generating OAM beams with topological charges of −2, −1, +1, and +2. Figure 3 shows the S11 amplitude phase polar coordinate statistical distribution obtained by simulating 12,000 sets of metasurface unit structures generated from random data at 1.5 THz and 1.6 THz. It can be seen that the reflection amplitude of the vast majority of metasurface units tends to approach 1, and the reflection phase is significantly concentrated in the range of 315°–45° (i.e. 0° ± 45°). This distribution truly reflects the inherent statistical manifold characteristics of the metasurface response space under unbiased random sampling, verifying the sample representativeness of the dataset in the low loss working range, and providing a benchmark distribution that conforms to statistical rules for subsequent neural network learning of physical laws.

Process of pixelization assisted design metasurfaces.
Figure 2. Process of pixelization assisted design metasurfaces.
The amplitude phase polar coordinate statistical distribution of the reflection coefficient S<sub>11</sub> of randomly generated pixelated metasurface units at **(a)** 1.5 THz, **(b)** 1.6 THz.
Figure 3. The amplitude phase polar coordinate statistical distribution of the reflection coefficient S11 of randomly generated pixelated metasurface units at (a) 1.5 THz, (b) 1.6 THz.

In order to select a suitable forward prediction network, a comparison was made between the residual convolutional network (ResCNN) and five architectures including CNN, U-Net, Transformer, and CNN Transformer on the dataset. The training loss and validation loss curves are shown in Figure 4. We can find that ResCNN can quickly decrease and stabilize compared to conventional CNN, U-Net, and CNN Transformer in the training process. However, the key to demonstrating the model’s generalization ability lies in the performance of the validation set: ResCNN’s validation loss converges to 0.015, significantly better than conventional CNN’s validation loss of 0.067, U-Net’s validation loss of 0.05, CNN Transformer’s validation loss of 0.054, and Transformer’s validation loss of 0.145. Based on the above dataset, a forward prediction neural network is constructed using a deep residual convolutional network (as shown in Figure 5a), with a binary image matrix of 16 × 16 elements on the metasurface as input parameters, to predict its corresponding S parameters (i.e., the real and imaginary parts of the complex reflection coefficient S11, a total of 201 frequency points, and a total of 402 dimensional outputs). After a layer of 3 × 3 convolution (i.e., step 1 and padding 1), the single channel geometric information is extended to 64 feature channels while maintaining the same spatial resolution. The feature extraction stage adopts a multi-level residual block stacking strategy (ResConvBlock), where each residual block contains two 3 × 3 convolutional layers, batch normalization, and ReLU activation functions, which alleviate the gradient vanishing problem through skip connections. The number of channels is gradually doubled from 64 to 128 to 256, abandoning traditional downsampling operations, and the feature map space size is always maintained at 16 × 16. At the end of feature extraction, adaptive global average pooling is used to aggregate multidimensional spatial features into one-dimensional vectors, and finally a 402-dimensional S-parameter prediction value is output through a linear layer.

Comparison of forward prediction networks, **(a)** training loss and **(b)** validation loss.
Figure 4. Comparison of forward prediction networks, (a) training loss and (b) validation loss.
The detailed structure of the proposed network, **(a)** Forward prediction neural network (deep residual convolutional network), **(b)** Generative inverse network.
Figure 5. The detailed structure of the proposed network, (a) Forward prediction neural network (deep residual convolutional network), (b) Generative inverse network.

We design a serial architecture generative inverse network as shown in Figure 5b, which inputs the predicted geometric structure into a pre-trained forward network with frozen parameters to reconstruct the corresponding S-parameters of the metasurface. The input of the inverse network is the concatenated vector of the real and imaginary parts of the target electromagnetic response S11, with a dimension of 402. In the feature encoding stage, a fully connected feedforward network is used to expand the input dimension to 1024. After each layer, batch normalization, ReLU activation function, and Dropout operation with a dropout rate of 0.1 are sequentially applied to enhance the model’s generalization ability. In the decoding stage, the 1024-dimensional bottleneck vector obtained from encoding is reshaped into a 64 × 4 × 4 tensor, and upsampled step by step to 16 × 16 through two sets of transposed convolutions. Then, a 3 × 3 standard convolutional layer is added, and the output value is strictly constrained within the [0,1] interval by the Sigmoid activation function to approximate the binary physical distribution of the metasurface element structure. Under the serial training framework, gradient backpropagation is performed by calculating the physical loss between the predicted S parameters and the target S parameters, achieving inverse design without the need for real geometric labels. The physical loss function is consistent with the form used for training the forward network, and includes both the mean square error of the real/imaginary parts and the mean square error of the amplitude. Table 1 shows the serial neural network hyperparameters. By using the AdamW optimizer and gradient pruning strategy, the gradient norm is limited to within 1.0, ensuring that the updated step size remains in a stable range and allowing the loss curve to smoothly decrease from the initial high position. The training process is combined with the cosine annealing learning rate scheduling algorithm. When the validation loss does not improve for 150 consecutive rounds, an early stop mechanism is triggered to terminate the training. This serial framework fully utilizes the differentiable physical consistency of the forward electromagnetic simulator, significantly improving the generalization ability of the inverse model in situations with few samples and multiple solutions, and achieving efficient and physically reliable structure generation. The calculation is completed on a computer equipped with Intel i9-13900k (with a clock speed of 3.0 GHz), 64 GB of memory, and NVIDIA GeForce RTX 4090 graphics card.

TABLE 1. Serial neural network hyperparameters.

HyperparametersForward predictionGenerative inverse
OptimizerAdamW (β1 = 0.5, β2 = 0.999)AdamW (β1 = 0.5, β2 = 0.999)
Learning rate1 × 10−35 × 10−4
Batch size6464
LR SchedulerCosine annealing Tmax = 1000, Tmin = 1 × 10−6Cosine annealing Tmax = 1000, tmin = 1 × 10−6
Dropout rate0.050.1

Model iterative training is achieved by minimizing the loss function, which is a physical loss function as shown in Equation 1:where α = 0.3 is the weight coefficient of the amplitude loss term, R ^ i and I ^ i are the real and imaginary parts of the predicted S11 spectrum, respectively. The corresponding true values are R i, and I i, N = 201 is the number of frequency points. The predicted amplitude A ^ i and true amplitude A i are defined as Equation 2:where ε = 10−8 is a numerical stability constant. This loss function simultaneously constrains the complex numerical vector and its modulus, ensuring that the S11 generated by the inverse network is consistent with the target response in both phase and amplitude. We can evaluate the training results of neural networks by calculating the training and validation set losses. Figure 6a shows that the training and validation losses of the forward prediction neural network model rapidly decrease at the beginning of training, until the early stop mechanism is triggered after 459 rounds of training, and the training loss reaches 0.005 and the validation loss reaches 0.015. Figure 6d gives that after 724 rounds of training, the generative inverse network triggered the early stop mechanism, with a training loss of 0.01 and a validation loss of 0.026. We randomly select two sets of data from the testing sets of forward predictive neural network and generative inverse network. Figures 6b,c illustrate the error between the electromagnetic response predicted by the forward network and the simulated input value, while Figures 6e,f show that the predicted value of the generative inverse network has an error close to zero compared to the target value. The results further confirm the accuracy of the forward prediction neural network and the generative inverse network prediction. In order to further test the capabilities of the designed forward predictive neural network and generative inverse network models, the MSE was used to perform statistical analysis on the real and imaginary parts of the entire test set. The definition of MSE is calculated by Equation 3:where y i and y ^ i represent the true and predicted values of the real and imaginary parts of the S11 spectrum, respectively, and N is the total number of samples. Figures 7a,b depict the number of test samples within a specific prediction error range using histograms, where the blue dashed line represents the mean square error. The average MSE of the test sets for forward predictive neural networks and serial neural network are 0.0017 and 0.0131, respectively. The MSE statistics between the S11 obtained by reverse generating predicted metasurface units for full wave CST simulation value and the original S11 in the test set is shown in Figure 7c. The full wave simulation error of the reverse predicted structure is 0.0125. The difference between this result and the closed-loop test error of the serial neural network in Figure 7b is only 0.0006. It confirmed the reliability of the serial neural network. In addition, the reason why the error of serial neural networks is much higher than that of forward networks is due to the existence of one to many mappings in the electromagnetic design space.

L p h y s i c s = 1 − α · i N ∑ i = 1 N R ^ i − R i 2 + I ^ i − I i 2 + α · i N ∑ i = 1 N A ^ i − A i 2 (1)

A ^ = R ^ 2 + I ^ 2 + ε , A = R 2 + I 2 + ε (2)

MSE = 1 N ∑ i = 1 N y i − y ^ i 2 (3)

**(a)** Training and validation loss of forward prediction network model, **(b,c)** Comparison of predicted values and simulated values of real and imaginary parts of forward prediction network model, **(d)** Training and validation loss of generative inverse network model, **(e,f)** Comparison of predicted values and target values of real and imaginary parts of generative inverse network model.
Figure 6. (a) Training and validation loss of forward prediction network model, (b,c) Comparison of predicted values and simulated values of real and imaginary parts of forward prediction network model, (d) Training and validation loss of generative inverse network model, (e,f) Comparison of predicted values and target values of real and imaginary parts of generative inverse network model.
**(a)** MSE statistics of the forward prediction network test set, **(b)** MSE statistics of the serial neural network, **(c)** MSE statistics between the S<sub>11</sub> obtained by reverse generating predicted metasurface units for full wave CST simulation value and the original S<sub>11</sub> in the test set.
Figure 7. (a) MSE statistics of the forward prediction network test set, (b) MSE statistics of the serial neural network, (c) MSE statistics between the S11 obtained by reverse generating predicted metasurface units for full wave CST simulation value and the original S11 in the test set.

3 Result analysis

Here, we selected 1.5–1.6 THz as an example frequency band, designed 8 kinds of elements with different phase responses, defined 8 phase curves with a phase difference of 45° between them and an amplitude close to 1, converted these curves into real and imaginary parts input models. The top view of the designed the pixelated metasurface elements are shown in Table 2. Figure 8a shows that all elements within the target frequency band maintain amplitudes greater than 0.85, indicating high terahertz wave energy transmission efficiency. From Figure 8b, it can be seen that the designed 8 kinds of elements achieve comprehensive coverage of 360° phase range within the target frequency band, with a phase difference of about 45° between adjacent elements. In order to further validate the performance of the constructed pixelated metasurface model and promote its practical application, a metasurfaces with topological charge numbers of ±1 and ±2 OAM beams were arranged using the 8 kinds of metasurface elements. The metasurface array has 24 × 24 elements, as shown in Figure 9. To achieve effective spatial separation of multiple terahertz beams, a specific gradient guided phase distribution was designed at the center frequency of the target frequency band (−25°, +35° in the x-axis direction, ±30° in the y-axis direction) to generate the final composite phase encoding, so that the terahertz beam can be symmetrically directed towards the predetermined direction.

TABLE 2. Eight kinds of the predicted metasurface units.

Elements number01234567
Phase (1.55 THz)33.3°78.3°122.6°167.1°−148.1°−105.6°−60.1°−12.8°
StructurePixel art grid with a symmetrical geometric pattern in blue and white squares, featuring a central open area, alternating filled and empty sections, and mirrored horizontal and vertical design.Pixel art design consisting of a white grid with light blue squares arranged symmetrically to form an abstract geometric pattern that features a central square, corner accents, and interconnected shapes.Grid pattern illustration featuring two mirrored geometric designs in light blue on a white background, each forming an abstract symmetrical shape composed of small squares arranged within the grid.Grid illustration featuring blue-filled squares arranged in a symmetrical, abstract pattern with a white background. The blue shapes form a border-like design with a square, pixelated structure.Grid illustration featuring a symmetrical blue and white pixel pattern. Blue squares create concentric rectangles and curved shapes resembling a stylized geometric motif centered on a white grid background.Grid illustration showing a symmetrical pixel pattern with light blue squares forming geometric shapes, including an X motif at each corner and horizontal bars in the center. White squares fill the background.Pixel art design on a square grid features a symmetrical pattern of light blue and white squares, forming an abstract geometric shape resembling a simplified space invader or stylized emblem.Grid-based illustration showing a symmetrical blue geometric pattern with rectangles and squares on a white background, concentrated at the center and corners with evenly spaced grid lines for alignment.
The amplitude and phase of the metasurface units, **(a)** The predicted amplitude of the metasurface units, **(b)** The predicted phase of the metasurface units.
Figure 8. The amplitude and phase of the metasurface units, (a) The predicted amplitude of the metasurface units, (b) The predicted phase of the metasurface units.
OAM multimode multiplexing and beam deflection phase superposition operation, **(a)** OAM topology charge is *l* = +1, **(b)** *l* = −1, **(c)** *l* = −2, **(d)** *l* = +2, **(e)** Phase is offset by 25° along the negative x-axis direction, **(f)** Offset by 35° along the positive *x*-axis direction, **(g)** Offset by 30° along the positive y-axis direction, **(h)** Offset by 30° along the negative y-axis direction, **(i)** *l* = ± 1 multiplexing distribution, **(j)** *l* = ±2 multiplexing distribution, **(k)** Superimposed phase distribution.
Figure 9. OAM multimode multiplexing and beam deflection phase superposition operation, (a) OAM topology charge is l = +1, (b) l = −1, (c) l = −2, (d) l = +2, (e) Phase is offset by 25° along the negative x-axis direction, (f) Offset by 35° along the positive x-axis direction, (g) Offset by 30° along the positive y-axis direction, (h) Offset by 30° along the negative y-axis direction, (i) l = ± 1 multiplexing distribution, (j) l = ±2 multiplexing distribution, (k) Superimposed phase distribution.

When the y-polarized plane wave is vertically incident, the amplitude of the far-field reflected electric field on the metasurface is shown in Figure 10. In the frequency range of 1.5 THz to 1.6THz, the metasurface generates OAM with topological charge numbers of +2, −2, +1, and −1. The OAM beam with a topological charge numbers of+2 is located at (φ = 270°,θ = 30°), the OAM beam with a topological charge numbers of −2 is located at (φ = 90°,θ = 30°), the OAM beam with a topological charge numbers of +1 is located at (φ = 180°,θ = 25°), and the OAM beam with a topological charge numbers of −1 is located at (φ = 0°,θ = 35°). The four OAM beams are basically consistent with the expected OAM mode, with a typical “donut” shaped electric field intensity radiation, as shown in Figure 11. The phase distribution of different OAM modes varies along concentric circles in the frequency range of 1.5 THz∼1.6 THz. For topological charge numbers of +1 and +2, the OAM phase distribution rotates counterclockwise and increases. For topological charge numbers of −1 and −2, the OAM phase distribution rotates clockwise and increases. The phase distribution of ±1 topological charge changes by 2π during one rotation of the electric field, while the phase distribution of ±2 topological charge changes by 4π during one rotation of the electric field. Figures 12, 13 show that when the metasurface array is expanded from 24 × 24 to 32 × 32, the phase singularity splitting phenomenon originally observed in the l = +2 completely disappears, and the far-field OAM ellipticity also decreases, restoring the single center singularity structure. When the quantization level changes to 4 bits (The 16 predicted metasurface units are shown in Table 3.), Figures 14, 15 undergo the aforementioned changes. The previously observed phase singularity splitting is caused by diffraction truncation effects resulting from both limited aperture and insufficient phase quantization accuracy. For the l = +2, the spiral phase needs to accumulate a phase change of 4π in azimuth, and the phase gradient is steeper. The wavefront structure is extremely sensitive to aperture size and phase discretization errors. Expanding the array to 32 × 32 increased the effective aperture and directly alleviated the truncation effect. Raising the quantization level to 4 bits can make the discrete phase closer to the ideal continuous distribution, effectively suppressing wavefront distortion caused by truncation. Both approaches enable the correct 4π spiral phase structure encoded to be fully presented, and the single central phase singularity is also restored accordingly. Figure 16 shows the mode purity of vortex beams at 1.55 THz for a 24 × 24 array with 3-bit metasurface under different topological charge numbers. The mode purities for topological charges l = −2, l = −1, l = 1, and l = 2 are 75.7%, 82.6%, 91.8%, and 78.6%, respectively. Figure 17 displays the mode purity of vortex beams at 1.55 THz for a 32 × 32 array with 3-bit metasurface under varying topological charges, with purities of 84.3%, 83.6%, 90.2%, and 82.8% for l = −2, l = −1, l = 1, and l = 2, respectively. Figure 18 shows the comparison of vortex beam mode purity between 24 × 24 and 32 × 32 arrays of 3-bit metasurface. Figure 19 presents the mode purity of vortex beams at 1.55 THz for a 24 × 24 array with 4-bit metasurface under different topological charges, showing purities of 77.3%, 82.7%, 83.8%, and 83.1% for l = −2, l = −1, l = 1, and l = 2, respectively. These results not only verify the OAM characteristics of the generated beam from both the amplitude and phase dimensions, but also indicate that the proposed pixelated metasurface and deep learning inverse design method can accurately generate complex wavefront devices that can simultaneously carry multiple propagation directions and multiple topological loads according to the target response, fully reflecting the strong potential of deep learning driven inverse design in multi-functional terahertz wavefront regulation. The methods used in this article and various representative deep learning metasurface design strategies are shown in Table 4. It can be seen that our work has certain advantages in training epochs and datasets when achieving the predicted loss value. This article involves terahertz OAM beam generation and mode purity detection. The difficulty of experimental verification lies not only in sample preparation, but also in accurate measurement of OAM mode purity. Mode purity testing requires mode decomposition of the complete complex amplitude information of terahertz waves. However, traditional terahertz wave intensity detection cannot distinguish different OAM modes with similar intensity distributions, and phase interferometry measurement is extremely sensitive to environmental and manufacturing errors, making the experiment difficult. We currently do not have the relevant experimental testing equipment and are unable to conduct related experimental testing.

Simulation amplitude of far-field electric field with orbital angular momentum (OAM) for a 24 × 24 array, **(a)** Far-field distribution at 1.5THz, **(b)** Far-field distribution at 1.55THz, **(c)** Far-field distribution at 1.6 THz.
Figure 10. Simulation amplitude of far-field electric field with orbital angular momentum (OAM) for a 24 × 24 array, (a) Far-field distribution at 1.5THz, (b) Far-field distribution at 1.55THz, (c) Far-field distribution at 1.6 THz.
Far field phase characteristics of orbital angular momentum (OAM) for a 24 × 24 array under 4-bit quantization at 1.5 THz∼1.6 THz.
Figure 11. Far field phase characteristics of orbital angular momentum (OAM) for a 24 × 24 array under 4-bit quantization at 1.5 THz∼1.6 THz.
Simulation amplitude of far-field electric field with orbital angular momentum (OAM) for a 32 × 32 array, **(a)** Far-field distribution at 1.5THz, **(b)** Far-field distribution at 1.55THz, **(c)** Far-field distribution at 1.6 THz.
Figure 12. Simulation amplitude of far-field electric field with orbital angular momentum (OAM) for a 32 × 32 array, (a) Far-field distribution at 1.5THz, (b) Far-field distribution at 1.55THz, (c) Far-field distribution at 1.6 THz.
Far field phase characteristics of orbital angular momentum (OAM) for a 32 × 32 array at 1.5 THz∼1.6 THz.
Figure 13. Far field phase characteristics of orbital angular momentum (OAM) for a 32 × 32 array at 1.5 THz∼1.6 THz.

TABLE 3. 16 predicted metasurface units.

Elements number01234567
Phase (1.55 THz)8.2°32.1°56.2°77.8°98.7°120.5°141.6°163.3°
StructurePixel art pattern features a symmetrical geometric design in light blue on a white grid, forming four curved motifs arranged around a central rectangular shape with distinct, evenly spaced squares.Grid-based illustration featuring a symmetric pattern of blue squares on a white background, creating a geometric design with repeating horizontal and vertical shapes across a 12 by 12 cell grid.Pixel art illustration featuring a symmetrical blue pattern on a white grid, reminiscent of an abstract insect or geometric shape, centered within a square layout composed of individual squares.Pixel art graphic features a blue geometric figure resembling an 8-bit space invader icon, placed symmetrically on a white grid background with clear square outlines.Pixel art design on a grid with light blue squares forming a symmetrical abstract pattern, featuring an X shape, two vertical central bars, and additional dots at the corners and mid-sides.Grid illustration showing a symmetrical geometric pattern with light blue squares on a white background, forming cross-like shapes and mirrored angles, centered within a square grid structure.Pixel art grid with blue squares forming a symmetrical pattern, featuring a central open square surrounded by angular shapes, resembling a stylized cross or abstract geometric design.Grid pattern illustration with light blue squares arranged symmetrically on a white background, forming a diamond shape with additional blue blocks along the top and bottom edges.
Elements number89101112131415
Phase (1.55 THz)174.6°−166.1°−144.1°−118.7°−97.9°−69.7°−44.4°−16.7°
StructureGrid-based design featuring a symmetrical blue geometric pattern with mirrored elements in all four quadrants, forming an abstract, decorative motif on a white background. Each square in the grid is clearly defined.Pixel art grid featuring a symmetrical blue geometric pattern resembling an abstract shape, centered on a white background with equally spaced squares, suggestive of a cross-stitch or embroidery template.Pixel art illustration on a grid depicts a stylized crab or insect with a rectangular body in the center, two claws extended to each side, and symmetrical patterns above and below the body in blue squares.Pixel art grid depicting a symmetrical butterfly or insect-like figure using blue squares on a white background, with patterned wings and body centered within a bordered square layout.Pixel art illustration featuring a symmetrical blue and white abstract pattern within a twenty by twenty grid, with blue squares forming geometric shapes against a white background.Symmetrical pixel art pattern on a 15 by 15 square grid showing a mirrored abstract design formed by blue squares on a white background, with shapes in each quadrant converging toward the center.Pixel art graphic showing a symmetrical circular pattern in light blue centered on a white grid. The design resembles a stylized mandala or geometric snowflake with repeating angular shapes.Pixel art pattern on a white grid with light blue squares arranged in a symmetrical geometric design, featuring mirrored shapes in each quadrant and a central white square.
Simulation amplitude of far-field electric field with orbital angular momentum (OAM) for a 24 × 24 array with 4-bit metasurface, **(a)** Far-field distribution at 1.5THz, **(b)** Far-field distribution at 1.55THz, **(c)** Far-field distribution at 1.6 THz.
Figure 14. Simulation amplitude of far-field electric field with orbital angular momentum (OAM) for a 24 × 24 array with 4-bit metasurface, (a) Far-field distribution at 1.5THz, (b) Far-field distribution at 1.55THz, (c) Far-field distribution at 1.6 THz.
Far field phase characteristics of orbital angular momentum (OAM) for a 24 × 24 array under 4-bit quantization at 1.5THz~1.6THz.
Figure 15. Far field phase characteristics of orbital angular momentum (OAM) for a 24 × 24 array under 4-bit quantization at 1.5THz~1.6THz.
Mode purity of vortex beams at 1.55 THz for a 24 × 24 array with 3-bit metasurface under different topological charges. **(a)** l = -2, **(b)** l = -1, **(c)** l = 1, **(d)** l = 2.
Figure 16. Mode purity of vortex beams at 1.55 THz for a 24 × 24 array with 3-bit metasurface under different topological charges. (a) l = -2, (b) l = -1, (c) l = 1, (d) l = 2.
Mode purity of vortex beams at 1.55 THz for a 32 × 32 array with 3-bit metasurface under different topological charges. **(a)** l = -2, **(b)** l = -1, **(c)** l = 1, **(d)** l = 2.
Figure 17. Mode purity of vortex beams at 1.55 THz for a 32 × 32 array with 3-bit metasurface under different topological charges. (a) l = -2, (b) l = -1, (c) l = 1, (d) l = 2.
Comparison of vortex beam mode purity between 24 × 24 and 32 × 32 arrays of 3-bit metasurface.
Figure 18. Comparison of vortex beam mode purity between 24 × 24 and 32 × 32 arrays of 3-bit metasurface.
Mode purity of vortex beams at 1.55 THz for a 24 × 24 array with 4-bit metasurface under different topological charges. **(a)** l = -2, **(b)** l = -1, **(c)** l = 1, **(d)** l = 2.
Figure 19. Mode purity of vortex beams at 1.55 THz for a 24 × 24 array with 4-bit metasurface under different topological charges. (a) l = -2, (b) l = -1, (c) l = 1, (d) l = 2.

TABLE 4. Comparison of this work with previously Reported studies.

RefsDesign methodNumber of tunable parametersTraining epochsLoss valueDataset
[19]Tandem neural network2300010–28281
[20]CNN64200010–123,000
[21]GAN5025003.1 × 10−320,000
This workTandem neural network49<6001.31 × 10−212,000

4 Conclusion

We used forward prediction residual convolutional networks and generative inverse metasurface network models to achieve forward prediction and inverse design of metasurfaces. The average MSE of the forward network test set is 0.0017, and the average MSE of the inverse network test set is 0.0131. It can achieve the design of the target electromagnetic response metasurface element structure, and realize the metasurface design with topological charge numbers of ±1 and ±2 OAM beams. Finally, the full wave simulation results indicate the effectiveness of the inverse design of high freedom pixelated metasurfaces.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Jun Zhu, Guangxi Normal University, China

Reviewed by: Eric Palmerduca, Princeton University, United States

Guirong Feng, Fuzhou University, China

Yangyang Zhou, Shantou University, China

Statements

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work is supported by the National Natural Science Foundation of China (62271460, 62435017), Zhejiang Provincial Natural Science Foundation of China under Grant No. LZ24F050005, and Fundamental Research Funds for the Provincial Universities of Zhejiang (2025YW108).

Data availability

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

RH: Formal Analysis, Software, Writing – original draft. JL: Conceptualization, Investigation, Supervision, Writing – review and editing. G-HQ: Formal Analysis, Writing – original draft.

References

  1. 1. Chen W, Wang Y, Cui Z, Sun G, Zhang K. Quasi-bound states in the continuum and sensing performance with polarization-selective in all-dielectric THz metasurface. IEEE Sens J (2024) 24: 6085 – 91. doi:10.1109/jsen.2024.3350039
  2. 2. Zhang C, Chen Z, Yi Z, Cheng S, Ahmad S, Tang C Broadband long-wave infrared metamaterial based on square-patterned Cr-based multilayer absorber. Mod Phys Lett B (2026) 40: 2650170. doi:10.1142/s0217984926501708
  3. 3. Luo H, Chen J, Liu L, Zhang Y, Yi Z, Cheng S Design and analysis of a compact quad-band polarization-insensitive metamaterial absorber based on coplanar nested resonators. Solid State Commun (2026) 416: 116502. doi:10.1016/j.ssc.2026.116502
  4. 4. Tuo C, Chen J, Liu L, Zhang Y, Yi Z, Cheng S Ultra-broadband mid-infrared metamaterial absorber based on cross-ring coupled cavities. Phys Status Solidi B (2026) 263: e70267. doi:10.1002/pssb.70267
  5. 5. Chen Z, Liu M, Cheng S, Wang J, Yi Y, Li B Bilayer graphene metasurface with dynamically reconfigurable terahertz perfect absorption. Curr Appl Phys (2025) 80: 282 – 90. doi:10.1016/j.cap.2025.10.001
  6. 6. Chen L, Cai Z, Yu X, Zhou Z, Dai Y, Guo J Full-space wave manipulations for multifunctional integration based on mechanically reconfigurable and stacked metasurface. Metamaterials (2026) 2: 17. doi:10.1038/s44455-026-00025-w
  7. 7. Gao Y, Chen W, Li F, Zhuang M, Yan Y, Wang J Meta‐attention deep learning for smart development of metasurface sensors. Adv Sci (2024) 11: 2405750. doi:10.1002/advs.202405750
  8. 8. Huang Y, Liu X, Liu M, Chen J, Du W, Liu Z. Deep learning-based inverse design of multi-functional metasurface absorbers. Opt Lett (2024) 49: 2733 – 6. doi:10.1364/OL.518786
  9. 9. Wang S, Chen G, Lu L, Zhu L. Inverse design and spectral reconstruction of computational multispectral metasurfaces using deep learning and genetic algorithms. Opt Express (2025) 33: 51485 – 98. doi:10.1364/OE.576847
  10. 10. Niu R, Si L, Zhu W, Bao X, Zhu K, Sun H Hybrid neural network-driven inverse design of metasurface for ultra-broadband RCS reduction. IEEE Trans Antennas Propag (2026) 74: 837 – 847. doi:10.1109/TAP.2025.3621515
  11. 11. Qi X, Xiao L, Lv H, Liu Y, Shao W. A machine learning-based inverse design method for frequency-selective surface microwave absorbers. Eng Appl Artif Intell (2025) 160: 111842. doi:10.1016/j.engappai.2025.111842
  12. 12. Ding Z, Su W, Ye L, Li W, Zhou Y, Tang B Deep learning based inverse design of metasurface absorber for maximizing solar spectral absorption. Sol Energy (2024) 271: 112449. doi:10.1016/j.solener.2024.112449
  13. 13. Yan Y, Wu Y, Wang Y, Li J, Hu J, Wang X. Deep-learning-based polarization-dependent switching metasurface in dual-band for optical communication. Nanophotonics (2025) 14: 5141 – 52. doi:10.1515/nanoph-2025-0370
  14. 14. Wu Q, Fan WH, Qin C, Jiang XQ. Inverse design of reconfigurable terahertz metasurface with wide phase modulation range via a high–precision CNN–Transformer framework. J Phys D: Appl Phys (2025) 58: 475103. doi:10.1088/1361-6463/ae1e89
  15. 15. Chen S, Hou Z, Wang J, Qiu Y, Li Z, Shen J Chiral metasurface vortex focusing in terahertz band based on deep learning. IEEE Photon Technol. Lett. (2023) 35: 637 – 40. doi:10.1109/lpt.2023.3264472
  16. 16. Guo Q, Hao F, Qu M, Su J, Li Z. Multiband multifunctional polarization converter based on reconfigurable metasurface. IEEE Antennas Wirel Propag Lett (2024) 23: 1241 – 5. doi:10.1109/lawp.2024.3350034
  17. 17. Yang H, Zheng S, Zhang H, He T, Li N, Yang Z Metasurface-based high-speed photonic THz OAM communication system. J Lightwave Technol (2024) 42: 5080 – 7. doi:10.1109/jlt.2024.3384477
  18. 18. Zhan J, Zhang Y, Li Q, Song Z. Broadband multifunctional metasurface for dynamic wavefront modulation. Opt Lasers Eng (2025) 194: 109193. doi:10.1016/j.optlaseng.2025.109193
  19. 19. Wei H, He X, Cao W. Spin-multiplexed metasurface inverse-design based on a bi-directional deep neural network for terahertz wavefront control. Optica (2025) 12: 505 – 17. doi:10.1364/optica.549420
  20. 20. Cheng J, Li R, Wang Y, Yuan Y, Wang X, Chang S. Inverse design of generic metasurfaces for multifunctional wavefront shaping based on deep neural networks. Opt Laser Technology (2023) 159: 109038. doi:10.1016/j.optlastec.2022.109038
  21. 21. Xia H, Chen S-L, Wang Y, Zhao Y, Jia H, Yang R Deep-learning-assisted intelligent design of terahertz hybrid-functional metasurfaces with freeform patterns. Opt Laser Technology (2025) 181: 112041. doi:10.1016/j.optlastec.2024.112041

Citation

Huang, R., et al. (2026). Inverse design of multi-mode OAM-generating pixelated THz metasurfaces. Frontiers in Physics, Volume 14. https://doi.org/
Received
Accepted
Published
Volume
14 - 2026

Copyright

Copyright © 2026 Huang, Li and Qiu.

License: https://creativecommons.org/licenses/by/4.0/

*Correspondence: Jiusheng Li

Article metrics

Views
0
Downloads
0
Citations
0
Shares
0

Related articles

Original ResearchAdvanced visualization workflows for gyrokinetics simulation dataRui Costa, William Hornsby, Daniel Kennedy, Lorenzo Zanisi, Nitesh Bhatia, Stanislas J. P. Pamela, Alejandra N. Gonzalez-Beltran, Robert AkersOriginal ResearchPatient-specific hemodynamic evaluation of fenestrated physician-modified endografts in complex abdominal aortic aneurysmsXiaolong Du, Shichen Jiang, Hongyue Zhang, Huiying Yu, Li Zhang, Kai Ren, Feng RanOriginal ResearchAgate-like stripe pattern formation via cyclic-voltage-driven reaction-transport -reaction processHisashi Hayashi, Kaho Yamada