
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)01820-0
10.1016/j.isci.2024.110595
110595
Article
Physics-driven unsupervised deep learning network for programmable metasurface-based beamforming
Bao Jianghan 123
Li Weihan 123
Huang Siqi 12
Yu Wen Ming 12
Liu Che cheliu@seu.edu.cn
12∗
Cui Tie Jun tjcui@seu.edu.cn
124∗∗
1 The State Key Laboratory of Millimeter Waves, Southeast University, Nanjing 210096, China
2 Institute of Electromagnetic Space, Southeast University, Nanjing 210096, China
∗ Corresponding author cheliu@seu.edu.cn
∗∗ Corresponding author tjcui@seu.edu.cn
3 These authors contributed equally

4 Lead contact

26 7 2024
20 9 2024
26 7 2024
27 9 11059517 5 2024
1 7 2024
24 7 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Programmable metasurfaces have garnered significant attention for their capacity to dynamically manipulate electromagnetic (EM) waves. In particular, the programmable metasurfaces offer to generate a wide range of EM beams when the appropriate digital coding patterns are designed. Traditionally, optimizing the coding patterns involves time-consuming nonlinear optimization algorithms due to the high computational complexity. In this study, we propose a physics-assisted deep learning (DL) model that can calculate the coding pattern in milliseconds, requiring only a simple depiction of the desired beam. An extended version of the macroscopic model for digital coding metasurface is introduced as the physics-driven component, which can compute the radiation pattern rapidly based on the provided coding pattern. The integration of the macroscopic model ensures to generate the physics-compliant coding designs. We validate the proposed method experimentally by measuring several coding patterns for both single-beam and dual-beam scenarios, which demonstrate good performance of beamforming.

Graphical abstract

Highlights

• A quick and accurate method for the far field of large-scale metasurface is proposed

• The physics-driven part eases significantly the burden of data preparation

• The model can give coding patterns for both single- and dual-beam cases in milliseconds

Applied sciences; Machine learning

Subject areas

Applied sciences
Machine learning
Published: July 26, 2024
==== Body
pmcIntroduction

Electromagnetic (EM) metamaterials have captured significant attention over the past decades due to their remarkable capacities to manipulate the EM waves.1,2,3 These materials exhibit fascinating properties, including negative refraction,3 perfect lensing,4 invisibility cloaking,5,6,7,8 concentrators,9 illusion optics devices,10,11 and others.12 Metasurfaces, the 2-dimensional counterparts of metamaterials, consist of periodic or quasi-periodic sub-wavelength macroscopic units on a plane, with their properties determined by the arrangement of these meta-units.12,13 However, the traditional metasurfaces are constrained by their inability to adapt their EM characteristics to new requirements once they are fabricated, which significantly limits their applications. To address this challenge, a novel digital programmable metasurface (DPM) has been introduced.14 DPM employs digital encoding of the meta-units, such as using 0 and 1 to represent units with phase values of 0° and 180°. The states of these units are controlled by field-programmable gate array (FPGA) hardware, enabling dynamic and flexible manipulation of the EM waves. DPM has demonstrated its potentials in various fields, including manipulating the amplitude,15 polarization16 and orbital angular momentum17 of EM wave, holography,18,19,20 and quantum photonics.21

Traditional metasurface design methods primarily fall into two categories: back-projection (BP)22 and nonlinear optimization.23 BP considers the far-field radiation as the Fourier transform of the current distribution, allowing phase arrangements to be obtained via inverse Fourier transforms when provided with the desired radiation pattern. However, the BP method struggles to provide accurate coding designs for complex beam requirements. In contrast, nonlinear optimization approaches iteratively minimize objective functions to generate optimal metasurface coding patterns, making them suitable for more complex scenarios. Nonetheless, the time-intensive nature of these iterations hinders real-time responses to rapidly changing demands, posing a challenge for real-time and reliable coding design.

In recent years, deep learning (DL) models have achieved significant success across various domains, such as speech recognition, image recognition, and automatic translation.24 Some studies have already explored the application of DL to EM problems,25,26 furthermore on the problem of metasurface design.27,28,29,30,31,32 However, these methods suffer from a significant inconvenience: they all require extensive data preparation using computationally complex and time-consuming traditional methods. Additionally, as a result of multiplicity-solution property of the iterative optimization algorithm, it requires selection of data to ensure that when these data are used to train the supervised DL methods, the model can converge properly.18 Thus, for the coding pattern design task of DPM, we strive to develop an unsupervised DL model that does not rely on dataset from nonlinear optimization or full-wave simulation to guide the training process. Drawing inspiration from the previous work on holographic metasurface design using physics-assisted generative adversarial networks,18 here we aim to reduce dependence on traditional methods by introducing more EM mechanisms.

The proposed model comprises two main components. The first component is ResNet,33 a widely used DL framework in computer vision. Its input images aim to imitate the radiation pattern of the desired beam. The output of ResNet represents the phase arrangement of the metasurface. The second component is physics-driven, where the 2-dimensional radiation pattern is calculated using an extended macroscopic model (EMM)34 based on the ResNet’s output. Making use of full-simulation results, the newly proposed macroscopic model of metasurface can give much more accurate radiation pattern than traditional array theory approach.35 However, its original version still requires many simulation results for large-scale metasurface. To tackle this issue, we make an assumption that the radiation of distinct metaunits only differs by a phase difference term. This extended version needs only several results obtained by simulation to give far-field radiation accurately and feasibly. The ResNet and EMM correspond to the encoder and decoder in the autoencoders36 in DL, respectively. The major difference lies in that the parameters of the decoder are usually trainable, whereas in our model, we replace it with non-trainable EMM. The introduction of this EM mechanism imposes additional constraints on the training of model parameters, thereby giving clear meaning in physics to the output, i.e., the coding distribution on the metasurface. The parameters of ResNet are optimized to minimize the difference between the output radiation pattern and the target image. Once this difference becomes sufficiently small, the ResNet output is considered as the reasonable metasurface coding pattern for the required beams. We explore the performance of different loss functions and select the more suitable cosine embedding loss (CEL) instead of the commonly used mean square error (MSE) loss. A schematic diagram of the proposed model is given in Figure 1.Figure 1 Conceptual illustration of the proposed approach

To better illustrate the effectiveness of our DL model, we compare the coding patterns designed by the proposed model with the results of nonlinear optimization approach, obtained by a modified binary particle swarm optimization (MBPSO) algorithm. The far-field radiation pattern generated by our method usually exhibits lower sidelobe. From the perspective of speed, the proposed DL method converges in 10 epochs, taking approximately 17 min for training. After convergence, it can provide the coding configurations for any single-beam or dual-beam scenarios in just a few milliseconds. Its capacity of rapidly forming beams extends its applicability to real-time beam scanning, microwave imaging, intelligent sensing, and so on. In contrast, the binary particle swarm optimization (BPSO) method requires about 10 min of optimization for each requirement on our device. The validity of the coding patterns output by the proposed method is also verified by experimental measurement. We probe the electric field intensities and phases for both single- and dual-beam cases with different beam orientations. The results show that our coding patterns are capable of generating beams at the required directions. The main contributions can be summarized as the following three aspects.(1) The original macroscopic model for the programmable metasurface is extended to adapt to the large-scale situation. Simulations of only four coding patterns are sufficient to provide accurate and fast results for any coding arrangements by using the extended version in our case.

(2) We propose a physics-driven DL model, enabling the model training without the need for extensive dataset preparation. Our model can design the coding patterns for both single- and dual-beam requirements in quasi real-time and the effectiveness of these designs is verified by experimental measurements.

(3) The BPSO algorithm is improved and compared with the proposed model. It is found that our model not only outperforms the traditional methods in terms of computation speed but also in the design effectiveness.

Digital coding metasurface

The hardware of programmable digital coding metasurface used for beamforming is presented in Figure 2, which has been used for target tracking before.37 A detailed description of the meta-unit is given in the following. As depicted in Figure 2A, the 1-bit coding unit comprises two PIN diodes connected to a central patch via two metal bars along the x and y directions. The unit’s geometric parameters are: a = 25 mm, b = 11.5 mm, c1 = 6 mm, c2 = 2.8 mm, d = 1.5 mm, r = 5 mm and β = 120°. Two dielectric substrates are employed in the unit, both of which are F4B material with the dielectric constant ϵr = 3 and tangent loss of tanδ = 0.003. The upper substrate has a height of h1 = 3 mm, while the lower one has a height of h2 = 1 mm. Since both diodes along the x and y directions have OFF and ON states, the entire meta-unit exhibits four distinct states. In this study, the diode along the y direction remains invariant, resulting in the meta-unit functioning as a 1-bit coding element. We use 0 and 1 to represent the OFF and ON states of the x axis diode. We simulate the element response under different frequencies using the commercial software, CST Microwave Studio, and the simulated results are presented in Figure 2C. Around the central frequency of 5.8 GHz, the normalized reflected amplitudes of both states are larger than 0.85 and the phase difference between the 0 and 1 states is approximately 180°. The complete aperture of the metasurface used in this work consists of 324 elements arranged in an 18 × 18 grid, as shown in Figure 2B. More information about the hardware of the metasurface can be found in Figure S4.Figure 2 The structure and performance of the metasurface

(A) The structure of the meta-unit and the bottom view of the coding element.

(B) Photo of the fabricated metasurface.

(C) The variations of reflection amplitude and phase of meta-unit at states 0 and 1 with respect to the frequency.

Deep learning method

Extended macroscopic model

Recently, a feasible and accurate method call macroscopic model for computing the scattered field of metasurface has been proposed in a study by Shao at al.34 Based on a small number of simulated results, this approach avoids the significant computational burden of full-wave simulation,38,39 which makes it convenient to be integrated into a DL framework. Compared with those methods where the far-field of meta-unit is considered to be omnidirectional or cosine,35,40,41 the macroscopic method is capable of giving much more accurate radiation pattern by effectively utilizing the results of simulations or measurements. By introducing the current space, which consists of the currents on the passive structures and tunable devices of metasurface, the equivalent current on the aperture can be approximated as power series of the reflection coefficients si=±1 of positive-intrinsic-negative (PIN) diodes. The scattered field of the metasurface can be easily obtained by utilizing the integral operation of the current with the scalar or dyadic Green’s function as the kernel. Since the integral operation is linear, the electric field is also expressed by following approximation, which is composed of power terms of si,(Equation 1) E=E0+∑m,nM,NE1m,nsm,n+12∑m,nM,N∑m′≠m,n′≠nM,NE2m,n,m′,n′sm,nsm′,n′+⋯

where the superscript and subscript of E indicate the indexes of meta-unit the order of electric field pattern, respectively. An appropriate number of coding configurations and corresponding radiation are needed to solve these field terms.

In the original version of macroscopic model, the terms of every meta-unit are regarded as totally different. Thus, for a large-scale metasurface, many pairs of coding and radiation patterns are required, which causes much inconvenience. To alleviate the burden of computation, we make an assumption that the only disparity between the radiations of distinct meta-unit is the phase difference caused by their position.35 If we only consider the zeroth-order and first-order expansion, the EMM gives the electric field as follows:(Equation 2) E(uˆ)=E0(uˆ)+∑m,nM,NE1(uˆ)sm,ne−jk(|rm,n−rf|−rm,n·uˆ)

where rm,n and rf represent the location of mn-th meta-unit and feed horn, respectively, k is the wave vector, and uˆ indicates the direction of observation. In this paper, the feed horn illuminates the metasurface from a position 45 cm away at a 45° tilt angle. Four coding patterns and their radiated field computed by CST Microwave Studio are used to obtain the expansion terms by least square method. The four coding configurations and corresponding radiated field can be found in Figure S1. For these four coding patterns under consideration, Equation 2 can be converted to be following matrix form:(Equation 3) XEEMM=Y

where Y=[Ef1(uˆ),Ef2(uˆ),Ef3(uˆ),Ef4(uˆ)]T is the simulated far-field radiation of the four coding configurations. EEMM=[E0(uˆ),E1(uˆ)]T represents the expansion terms. X∈C4×2 and its i-th row is [1,∑m,nM,Nsi,m,ne−jk(|rm,n−rf|−rm,n·uˆ)], where si,m,n denotes the coding of the mn-th meta-unit in the i-th coding pattern. The least square method gives(Equation 4) EEMM=(XTX)−1XTY

The results of expansion terms are shown in Figure 3A. E0 indicates the part irrelevant to the change of coding pattern in electric field of the whole metasurface while E1 represents the radiation of just one meta-unit. E0 is significantly stronger than E1 and there are some strips in E0. With these terms, far-field radiation of any coding arrangement can be obtained nearly instantly. The comparison between the radiation patterns in uv-plane given by EMM and CST simulation of the same coding pattern, which is not included in the four coding configurations used to calculate the expansion terms, is given in Figure 3B. The high similarity between them suggests the effectiveness of EMM. It should also be emphasized that, the field in the dotted line rectangle is caused by specular reflection of metasurface, which is hardly effected by the change of coding pattern. The higher-order expansion terms and the poor coding design performance when we consider higher-order terms can be found in Figures S2 and S3, respectively. This indicates that the phase difference term only may not be sufficient to completely reveal the nature of mutual interactions at different positions.Figure 3 The results related to EMM

(A) The zeroth-order and first-order terms obtained by EMM. The first row is magnitude, while the second indicates phase. Eθ and Eϕ represent the electric field in the direction perpendicular to θ and ϕ direction, respectively.

(B) Comparison between the far-field radiation given by EMM and CST simulation.

Dataset

There have been many methods using DL approaches to design metasurfaces.27,28,29,30,31,32 However, a major drawback of these methods is the cumbersome preparation of datasets, which often involves full-wave simulations or nonlinear optimizations. For example, Shan et al.28 firstly utilizes genetic algorithm to optimize numerous coding patterns for the desired beamforming and then use them as targets when training the neural networks. These results given by the nonlinear optimizations also require selection to ensure the convergence of DL training. Because nonlinear optimization often involves a stochastic optimization process, it may optimize completely different coding patterns for 2 closely similar targets, which can easily lead to non-convergence when training the neural networks with such data.18 The aforementioned methods require the datasets to be produced by rigorous computations and selections because they make use of the supervised learning, which necessitates precise guidance during the training process of neural network. By employing an unsupervised learning strategy, our model does not demand datasets produced by time-consuming computation. Our target patterns do not have to perfectly match the actual field distribution. We only need to specify in which region of uv-plane we wish the beam to appear. In order to better match the shape of beam in uv-plane, 2-dimensional Gaussian distribution is used to mark the zone illuminated by the desired beam. The center of the distribution (cu, cv) is given by(Equation 5) cu=Lsinθ0cosϕ0/2,cv=Lsinθ0sinϕ0/2

where θ0 and ϕ0 are the angles that characterized the direction of the beam, and the origin of coordinates is set at the geometric center of pattern. L represents the size of target images here and is 1,024. σ of the Gaussian distribution is set to be 30. Illustration of aforementioned parameters and four example target patterns are given in Figure 4. After the Gaussian distribution is generated, all the pixels will be normalized to ensure that the value at the center is 1. For the dual beam scenario, it’s only necessary to superpose the target patterns of two single beams. During training, the possible elevation angle of required beam θ0 is set to range from 0° to 45°, and azimuth angle ϕ0 starts from 0° to 360°. Their steps are 3°. We use both single-beam and dual-beam target patterns to train the DL model, so that it can handle both scenarios simultaneously.Figure 4 Several examples of target patterns

The radiation patterns in this figure share the same colorbar.

Architecture of proposed model

A comprehensive depiction of the proposed model is presented in Figure 5. Leveraging the ResNet framework, which is one of the most successful approaches in computer vision,33 we aim to generate the metasurface coding based on the desired radiation pattern. Our implementation involves ResNet34, consisting of 16 residual blocks equipped with instance normalization layers42 and ReLU activation functions.43 The target radiation pattern is first reshaped into 64 × 64 images, which are subsequently fed into ResNet34. The final activation function applied in ResNet34 is composed of sigmoid function S and it is given as(Equation 6) S(x)=11+e−x

(Equation 7) output=2S(x)−1

Figure 5 The architecture of the DL model in this work

We use “Conv” to represent the 2-dimensional convolution layer, of which the kernel size is also given in the figure. “FC” indicates fully connected layer and 324 is its output dimension. “CEL” stands for cosine embedding loss. The white and green dashed boxes in target pattern and radiation pattern given by EMM mark the corresponding beams. The radiation patterns in this figure share the same colorbar.

Sigmoid activation ensures that the values of each output unit fall in the range from 0 to 1. Therefore, the output of ResNet34 is mapped into [−1, 1]. The output of ResNet34 is configured as a 1 × 324 vector and is then reshaped to form an 18 × 18 coding pattern. It should be noticed that, the phase distribution obtained by ResNet34 consists of continuous phase, different from the real 1-bit metasurface, which has only 2 possible phases for each meta-unit. However, the results show that this approximation only has minor impact on the effectiveness of the proposed model, and the detailed information will be described in the discussion section. The predicted radiation pattern is calculated from this coding pattern using Equation 2, which is grounded in the EM wave theory, rendering this approach physics-driven. The difference between the prediction and the target pattern serves as the basis for training the DL model once an appropriate loss function is selected. Here we choose the CEL rather than the commonly used MSE. The expression of the former is as follows:(Equation 8) CEL(x1,x2)=1−cos(x1,x2)

where x1 and x2 denote the flattened 1-dimensional vectors of the target and output patterns, respectively. The cosine similarity is obtained via inner production of 2 vectors. Our experiments show that the MSE loss function is nearly ineffective for our task, as it fails to optimize beams with different orientations. We will also elaborate on this in Section V-D.

Binary particle swarm optimization method

To illustrate the efficiency of our DL method in designing the metasurface coding pattern, we compare the coding patterns generated by the proposed approach and those given by the BPSO. Particle swarm optimization (PSO) is a nature-inspired evolutionary optimization algorithm44 that simulates the movement and collaboration of a group of individuals, referred as particles, in a multidimensional search space. The position of a particle corresponds to a solution in the search space, and its movement is determined by both its own experience and the collective behavior of the swarm. In its binary version, i.e., BPSO, the particles consist of binary codes 0 and 1, which corresponds to the phase of elements in the problem of metasurface design in this paper. In the initial form of BPSO method, the j-th element of vi, which is the speed of the i-th particle on iteration k + 1 are updated by the following equations:(Equation 9) vijk+1=wvijk+c1·rand1·(pbestijk−xijk)+c2·rand2·(gbestjk−xijk)

where w is the inertia coefficient, c1 and c2 are acceleration constants, pbestijk is the j-th element of the optimal position for the i-th particle and gbestjk is the j-th element of global best position on iteration k. The speed vij represents the probability of xij taking the value of 0 and 1. BPSO uses sigmoid function S (Equation 6) to transform the value of speed into [0, 1]. Each element of i-th particle on next iteration is given by(Equation 10) xijk+1={1rand<S(vijk)0rand≥S(vijk)

where rand represents a random number. The original version of BPSO often suffers trapping in local minima,45 so several modifications are introduced here to enhance its performance. The first one is variable c1 and c2 with respect to iteration k to adjust the relative strength of self- and social-learning.46 The large c1 and small c2 at the beginning of optimization would strengthen the learning ability of each particle. The small c1 and large c2 near the end of optimization would let the particles focus on social-learning to get to global optimal solution. c1 and c2 at every iteration k are given by(Equation 11) c1(k)=c1,ini+c1,fin−c1,inikmaxk

(Equation 12) c2(k)=c2,ini+c2,fin−c2,inikmaxk

where c1,ini and c2,ini are initial values of c1 and c2, while c1,fin and c2,fin are final values for them. kmax denotes the total iteration number. The second modification is the introduction of nonlinear weight for each particle wi,46 which is defined by(Equation 13) wi={wmin+(wmax−wmin)(fi−fmin)favg−fmin,fi≤favgwmax,fi>favg

where wmax and wmin denote the maximum and minimum of w, respectively. fi is the value of fitness function for the i-th particle. favg is the average value of all particles and fmin is the minimum among them. When the fitness values verge to the same, the increasing w will enhance the global searching ability.

The last modification of BPSO is the replacement of the Sigmoid function (Equation 6) with another transfer function. According to Mirjalili et al.,45 a V-shaped one is selected to achieve improved performance. The new transfer function is given by(Equation 14) T(vij)={1−21+e−λvij,vij≤021+e−λvij−1,vij>0

where λ is a parameter to adjust the shape of T. The update strategy of position for vij ≤ 0 is(Equation 15) xij={0,rand≤T(vij)xij,rand>T(vij)

For vij > 0, we have(Equation 16) xij={1,rand≤T(vij)xij,rand>T(vij)

In summary, our modified BPSO (MBPSO) method is a combination of aforementioned 3 modifications and the initial BPSO approach. In the following part of this paper, we will compare the radiation patterns corresponding to coding patterns given by our DL model and MBPSO, as well as the speed of coding generation of these 2 approaches.

Results

Results of DL and MBPSO methods

In this section, we firstly present some details in the training process of the neural network. The training data contain both cases of single beam and dual beams. For the latter scenario, the difference between the θ angles of two beams ranges from 0° to 42° and that of ϕ ranges from 21° to 93°. The steps are still 3°. The feed horn is located at the direction of θ = 45°, ϕ = 0° and the distance from the metasurface is 50 cm. We use Adam optimizer47 to update the trainable weights in the neural networks and the learning rate is set to be 8 × 10−4. The batch size is 32 and training the DL network for 10 epochs is sufficient to obtain a well-converged model. Our code is constructed on the Pytorch platform and runs on a graphics processing unit (GPU) (Nvidia RTX A5000).

Some coding patterns designed by our method for different beam requirements are shown in Figure 6. The target patterns used as input of proposed DL model are listed in the first column along with the values of θ and ϕ of desired beam orientations. For every beam direction setting, the coding pattern designed by our model and the corresponding far-field radiation in uv-plane computed by EMM are given in the green box in the same row. The orange squares in the coding configurations represent the meta-units in 1 state, while the white squares are those in 0 state. It is clearly shown that the coding patterns given by our approach have close radiation profiles to the target patterns for both single-beam and dual-beam situations with distinct beam directions.Figure 6 Several coding designs given by proposed DL model and MBPSO algorithm

All the radiation patterns are normalized and share the same colorbar. Dashed boxes represent the corresponding beams.

The MBPSO algorithm is implemented via Python 3.9 platform and runs on a central processing unit (CPU) (AMD Ryzen 9 5950X 16-Core Processor). The parameters to control the optimization process are set as follows: wmin = 0.5, wmax = 0.99, c1,ini = 0.1, c1,fin = 2.0, c2,ini = 0.5, c2,fin = 2.0, λ = 1. The initial state of the particles is given by BP algorithm22 to accelerate the optimization. In contrast, the weights in our DL model are randomly initialized. For better comparison of the effects of DL method and MBPSO, we optimize the same target patterns using these two distinct approaches, and the results of MBPSO are presented in the blue box in Figure 6. It can be seen that the far-field radiation generated by the coding given by MBPSO have a similar general behavior to the target patterns. However, the sidelobe is obviously stronger than the designs of proposed DL model. Moreover, the MBPSO method has to optimize each objective from scratch, and our code takes about 10 min to run the optimization of 2,000 steps for 500 particles to achieve satisfactory results, and there also exists possibility to get completely unusable coding designs. In contrast, once the DL model converges, it can provide results for any objective within a few milliseconds, which brings a significant improvement in efficiency. A training process of 10 epochs of our DL model only needs around 17 min on our GPU.

Experimental verification

In order to further demonstrate the effectiveness of our method, we experimentally measure the coding patterns provided by the physics-driven DL model under different beam requirements. The experimental setup of the measurement is established in a microwave anechoic chamber, as illustrated in Figure 7.Figure 7 Experimental setup to measure the beams

The measured electric fields of 4 different coding designs are illustrated in Figure 8. The first column is the coding configurations designed by our physics-driven DL model, the second column presents the target patterns, the third one gives the corresponding radiation patterns computed by EMM, the fourth and fifth columns are the measured strength and phase of electric field, respectively. The measured strength is normalized to value in [0, 1] and all the field intensity diagrams share the same colorbar. Due to the high cost of scanning the entire 2-dimensional uv-plane, the probe is set to scan a rectangle area where θ is about to be smaller than 25°. And the scanned zone is not exactly at the center of the radiation pattern as a result of the limitation of probe orbit. The scanned area is indicated by red dashed box in EMM radiation patterns.Figure 8 Electric field strength and normalized phase measurements of four coding designs

Corresponding beams in target patterns, EMM radiations, and measured magnitudes are marked by white and green dashed boxes. The radiation patterns in this figure show the same colorbar.

The four measured coding patterns include four typical beamforming scenarios. The first one is intended to generate a single beam pointing directly forward. The result shows that our design model is able to give coding to compensate the phase difference caused by a horn antenna with oblique incidence. In the second case, we attempt to generate another beam with slant orientation additionally to that in the first case. The difference between these two coding configurations shows that modifying the stripe shape of the coding is sufficient to generate a second beam without affecting the performance of the first beam. We attempt to address a more challenging problem with the third coding pattern, which is required to form two beams in oblique direction. The measured data indicate that the strongest field is located as demanded directions. We then increase the value of θ for one of the beams in the last scenario. In the desired orientations, we have strong field; however, the beam in the right is stronger than that in the left. The right beam seems to be strengthened by the specular reflection of metasurface, as the horn illuminates it from the left direction of the figure. Generally speaking, the measurement demonstrates that the proposed model is capable of designing reliable coding pattern for both single- and dual-beam requirements.

Example of beam scanning

To better illustrate the potential application of our model, we give an example of tracking the target under different situation. As shown in Figure 9, we attempt to sequentially change the coding pattern of metasurface to adjust the orientation of beams for tracing the toy cars at different positions. In Figure 9A, the schematic diagram of different single beam for single target is given. The corresponding coding patterns and measured electric field are given in Figure 9B. For a more complex two targets case, the schematic diagram is shown in Figure 9C and experimental measurement is in Figure 9D. Since our DL model can output reliable coding configurations in milliseconds, the aforementioned results show that the beams formed by our coding patterns are able to catch plural targets in different locations.Figure 9 Example of beam scanning to trace the cars

(A) Illustration of the scenario of tracking single target using single beam.

(B) The coding patterns to realize the beams in (A) and the corresponding measured electric fields.

(C) Illustration of the scenario of tracking two targets using dual beams.

(D) The coding patterns to realize the beams in (C) and the corresponding measured electric fields. The desired beams are indicated by dashed boxes. The radiation patterns in this figure share the same colorbar.

Discussion

It should be emphasized that, the phase configuration given directly by the DL model contains continuous value in [−1, 1]. The 1-bit coding patterns in Figures 6 and 8 are obtained by discretization of the output with continuous phase and the predicted radiation patterns are given by the binary coding. For a converged DL model, we compare the far-field radiation patterns of its output before and after phase discretization, and found that there exists minor difference between the two settings (Figure 10A). The discretization only slightly increases the sidelobes. We plot a bar chart of the numerical values of the continuous phases given by a converged DL model for randomly sampled 720 target patterns as in Figure 10B. It can be observed from the chart that the distribution of continuous phase values concentrates around −1 and 1, resembling a 1-bit discrete distribution. This gives an explanation to the minor impact of this approximation.Figure 10 Illustration for the discretization of DL model output and choice of loss function

(A) In the upper row, we present an example of the continuous phase distribution directly output by a converged DL model and the corresponding radiation pattern obtained by EMM. In the lower row, we give the coding pattern after discretization and its far-field radiation.

(B) Distribution of phase output by a converged DL model.

(C) Optimization process of proposed model using different kinds of loss function. For the two beams in both cases, we have θ1 = 0°, ϕ1 = 0° and θ2 = 27°, ϕ2 = 357°. The number on the top of every column is the number of batches that have been used to update the model weights. Corresponding beams are marked by green and white dashed boxes. The radiation patterns in (C) share the same colorbar.

Now, we attempt to illustrate the inefficiency of MSE. We choose the same target patterns, with one beam at the center and another with oblique orientation, for both CEL and MSE and observe their optimization process. The far-field radiations based on the output by our DL model at several sampled training steps are presented in Figure 10C. The numbers of batches indicate how many batches have been utilized to optimize the model parameters. It can be seen that, for this dual-beam case, both of the loss functions can optimize the general profile of two beams. However, the CEL exhibits much quicker convergence of training. After 100 batches of parameter optimizing, the CEL is able to give two strong and clear beams. In contrast, the beams given by MSE are still weak and unbalanced. This indicates that CEL is the more appropriate choice.

Conclusion

We propose a physics-assisted DL model for the beamforming task based on metasurface in this study. This method does not require any paired training samples obtained through traditional approaches and also eliminates the need for selecting data to ensure training convergence. The simple target patterns generated based on the information of demanded beams are sufficient to train our model. Compared with commonly used non-linear optimization method (MBPSO in this paper), our DL method is capable of providing quasi real-time coding design once the training process is converged. The experimental measurements also validate the reliability of the coding patterns output of our model for both single- and dual-beam cases. The ability to quickly obtain the coding for the desired beams enables our method to be used in real-time beam scanning, microwave imaging, intelligent sensing, and other fields.

Limitations of the study

The EMM proposed in this work does not contain the coupling effect between the meta-units. If we include the mutual interactions in the EMM, the performance of the coding pattern designed by our method would be poor (refer to Figures S2 and S3 for details). This suggests that the difference between the coupling effects at distinct positions cannot be fully represented by the phase difference term. Therefore, further work is needed to address the mutual interaction properly.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Software and algorithms	
	
CST Studio Suite 2022	Dassault Systèmes	https://www.3ds.com	
Pytorch	Meta AI	https://www.pytorch.org	

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Tie Jun Cui (tjcui@seu.edu.cn).

Materials availability

This study did not generate new unique reagents.

Data and code availability

• Data reported in this paper will be shared by the lead contact upon request.

• The Pytorch model of the proposed method is publicly available via https://github.com/BaoJianghan/Physics_Driven_DL_DPM.

• Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Experimental model and study participant details

The CST Microwave Studio software has been employed to analyze the far-field pattern of the metasurface. The DL model is implemented on the Pytorch platform.

Method details

Simulate the metasurface with the four coding patterns listed in supplemental information (or other coding configurations). Apply the Time Domain Solver. Set boundary conditions as the open at x- and y-directions and the expanded open at z-direction. Set port excitation at the input end of the horn. Get the expansion terms of far-field electric pattern with the proposed EMM method. If the matrix is singular when using the least square method, turn to the pseudo inverse. Construct the DL model on Pytorch as section III C suggests. Set learning rate as 8 × 10−4, batch size as 32, and train the network for 10 epochs.

Quantification and statistical analysis

The simulation data are produced by CST Microwave Studio software. Figures shown in the main text were produced by Python Matplotlib package.

Additional resources

Any additional information about the simulation and data reported in this paper is available from the lead contact on request.

Supplemental information

Document S1. Figures S1–S4

Acknowledgments

This work was supported by the 10.13039/501100001809 National Natural Science Foundation of China (grant no. 62288101 ) and the Fundamental Research Funds for the Central Universities under grant 2242023K5002 .

Author contributions

Investigation, J.B., W.L., and C.L.; software, J.B., S.H., and C.L.; formal analysis, J.B. and C.L; resources, W.M.Y, C.L., and T.J.C.; writing – original draft, J.B. and C.L; writing – review & editing, W.M.Y, C.L., and T.J.C.; supervision, C.L. and T.J.C.; project administration, C.L. and T.J.C.; funding acquisition, C.L. and T.J.C.

Declaration of interests

The authors declare no competing interests.

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110595.
==== Refs
References

1 Zheludev N.I. Kivshar Y.S. From metamaterials to metadevices Nat. Mater. 11 2012 917 924 10.1038/nmat3431 23089997
2 Pendry J.B. Schurig D. Smith D.R. Controlling electromagnetic fields Science 312 2006 1780 1782 10.1126/science.1125907 16728597
3 Yu N. Genevet P. Kats M.A. Aieta F. Tetienne J.-P. Capasso F. Gaburro Z. Light propagation with phase discontinuities: Generalized laws of reflection and refraction Science 334 2011 333 337 10.1126/science.1210713 21885733
4 Pendry J.B. Negative refraction makes a perfect lens Phys. Rev. Lett. 85 2000 3966 3969 10.1103/PhysRevLett.85.3966 11041972
5 Schurig D. Mock J.J. Justice B.J. Cummer S.A. Pendry J.B. Starr A.F. Smith D.R. Metamaterial electromagnetic cloak at microwave frequencies Science 314 2006 977 980 10.1126/science.1133628 17053110
6 Li J. Pendry J.B. Hiding under the carpet: A new strategy for cloaking Phys. Rev. Lett. 101 2008 203901 10.1103/PhysRevLett.101.203901
7 Liu R. Ji C. Mock J.J. Chin J.Y. Cui T.J. Smith D.R. Broadband ground-plane cloak Science 323 2009 366 369 10.1126/science.1166949 19150842
8 Ma Q. Mei Z.L. Zhu S.K. Jin T.Y. Cui T.J. Experiments on active cloaking and illusion for laplace equation Phys. Rev. Lett. 111 2013 173901 10.1103/PhysRevLett.111.173901
9 Jiang W.X. Cui T.J. Cheng Q. Chin J.Y. Yang X.M. Liu R. Smith D.R. Design of arbitrarily shaped concentrators based on conformally optical transformation of nonuniform rational B-spline surfaces Appl. Phys. Lett. 92 2008 264101 10.1063/1.2951485
10 Lai Y. Ng J. Chen H. Han D. Xiao J. Zhang Z.-Q. Chan C.T. Illusion optics: The optical transformation of an object into another object Phys. Rev. Lett. 102 2009 253902 10.1103/PhysRevLett.102.253902
11 Chen L. Ma Q. Nie Q.F. Hong Q.R. Cui H.Y. Ruan Y. Cui T.J. Dual-polarization programmable metasurface modulator for near-field information encoding and transmission Photon. Res. 9 2021 116 10.1364/PRJ.412052
12 Cui T.J. Smith D. Liu R. Metamaterials: Theory, Design, and Applications 2009 Springer
13 Li L. Shi Y. Cui T.J. Electromagnetic Metamaterials and Metasurfaces: From Theory to Applications 2024 Springer
14 Cui T.J. Qi M.Q. Wan X. Zhao J. Cheng Q. Coding metamaterials, digital metamaterials and programmable metamaterials Light Sci. Appl. 3 2014 e218 10.1038/lsa.2014.99
15 Wu R.Y. Zhang L. Bao L. Wu L.W. Ma Q. Bai G.D. Wu H.T. Cui T.J. Digital metasurface with phase code and reflection–transmission amplitude code for flexible full-space electromagnetic manipulations Adv. Opt. Mater. 7 2019 1801429 10.1002/adom.201801429
16 Ma Q. Hong Q.R. Bai G.D. Jing H.B. Wu R.Y. Bao L. Cheng Q. Cui T.J. Editing arbitrarily linear polarizations using programmable metasurface Phys. Rev. Appl. 13 2020 021003 10.1103/PhysRevApplied.13.021003
17 Ding G. Chen K. Luo X. Zhao J. Jiang T. Feng Y. Dual-helicity decoupled coding metasurface for independent spin-to-orbital angular momentum conversion Phys. Rev. Appl. 11 2019 044043 10.1103/PhysRevApplied.11.044043
18 Liu C. Yu W.M. Ma Q. Li L. Cui T.J. Intelligent coding metasurface holograms by physics-assisted unsupervised generative adversarial network Photon. Res. 9 2021 B159 10.1364/PRJ.416287
19 Huang L. Zhang S. Zentgraf T. Metasurface holography: from fundamentals to applications Nanophotonics 7 2018 1169 1190 10.1515/nanoph-2017-0118
20 Ni X. Kildishev A.V. Shalaev V.M. Metasurface holograms for visible light Nat. Commun. 4 2013 2807 10.1038/ncomms3807
21 Solntsev A.S. Agarwal G.S. Kivshar Y.S. Metasurfaces for quantum photonics Nat. Photonics 15 2021 327 336 10.1038/s41566-021-00793-z
22 Bhattacharyya A.K. Phased Array Antennas: Floquet Analysis, Synthesis, BFNs and Active Array Systems 2006 Wiley
23 Isernia T. Massa A. Morabito A.F. Rocca P. On the optimal synthesis of phase-only reconfigurable antenna arrays Proceedings of the 5th European Conference on Antennas and Propagation 2011 EUCAP 2074 2077
24 LeCun Y. Bengio Y. Hinton G. Deep learning Nature 521 2015 436 444 10.1038/nature14539 26017442
25 Shan T. Tang W. Dang X. Li M. Yang F. Xu S. Wu J. Study on a fast solver for poisson’s equation based on deep learning technique IEEE Trans. Antennas Propag. 68 2020 6725 6733 10.1109/TAP.2020.2985172
26 Liu C. Zhang H. Li L. Cui T.J. Towards intelligent electromagnetic inverse scattering using deep learning techniques and information metasurfaces IEEE J. Microw. 3 2023 509 522 10.1109/JMW.2022.3225999
27 Qiu T. Shi X. Wang J. Li Y. Qu S. Cheng Q. Cui T. Sui S. Deep learning: A rapid and efficient route to automatic metasurface design Adv. Sci. 6 2019 1900128 10.1002/advs.201900128
28 Shan T. Pan X. Li M. Xu S. Yang F. Coding programmable metasurfaces based on deep learning techniques IEEE J. Emerg. Sel. Top. Circuits Syst. 10 2020 114 125 10.1109/JETCAS.2020.2972764
29 Niu C. Phaneuf M. Qiu T. Mojabi P. A deep learning-based approach to design metasurfaces from desired far-field specifications IEEE Open J. Antennas Propag. 4 2023 641 653 10.1109/OJAP.2023.3292108
30 Fu J. Zhang Y. Dou Z. Yang Z. Liu M. Zhang H. Rapid deep-learning-assisted design method for 2-bit coding metasurfaces Appl. Opt. 62 2023 3502 3511 10.1364/AO.487867 37132852
31 Qu K. Chen K. Hu Q. Zhao J. Jiang T. Feng Y. Deep-learning-assisted inverse design of dualspin/frequency metasurface for quad-channel off-axis vortices multiplexing Adv. Photon. Nexus 2 2023 016010 10.1117/1.APN.2.1.016010
32 Kiani M. Kiani J. Zolfaghari M. Conditional generative adversarial networks for inverse design of multifunctional metasurfaces Advanced Photonics Research 3 2022 2200110 10.1002/adpr.202200110
33 He K. Zhang X. Ren S. Sun J. Deep residual learning for image recognition IEEE Conference on Computer Vision and Pattern Recognition (CVPR) 2016 770 10.1109/CVPR.2016.90
34 Shao R.W. Wu J.W. Wang Z.X. Xu H. Yang H.Q. Cheng Q. Cui T.J. Macroscopic model and statistical model to characterize electromagnetic information of a digital coding metasurface Natl. Sci. Rev. 11 2024 nwad299 10.1093/nsr/nwad299
35 Nayeri P. Yang F. Elsherbeni A.Z. Radiation analysis techniques Reflectarray Antennas: Theory, Designs, and Applications Chap. 4 2018 John Wiley & Sons, Ltd 79 111
36 Bank D. Koenigstein N. Giryes R. Autoencoders Preprint at: arXiv 2020 10.48550/arXiv.2003.05991
37 Li W. Ma Q. Liu C. Zhang Y. Wu X. Wang J. Gao S. Qiu T. Liu T. Xiao Q. Intelligent metasurface system for automatic tracking of moving targets and wireless communications based on computer vision Nat. Commun. 14 2023 989 10.1038/s41467-023-36645-3 36813789
38 Francavilla M.A. Martini E. Maci S. Vecchi G. On the numerical simulation of metasurfaces with impedance boundary condition integral equations IEEE Trans. Antennas Propag. 63 2015 2153 2161 10.1109/TAP.2015.2407372
39 Bodehou M. González-Ovejero D. Craeye C. Huynen I. Method of moments simulation of modulated metasurface antennas with a set of orthogonal entiredomain basis functions IEEE Trans. Antennas Propag. 67 2019 1119 1130 10.1109/TAP.2018.2880075
40 Najafi M. Jamali V. Schober R. Poor H.V. Physics-based modeling and scalable optimization of large intelligent reflecting surfaces IEEE Trans. Commun. 69 2021 2673 2691 10.1109/TCOMM.2020.3047098
41 Di Renzo M. Danufane F.H. Tretyakov S. Communication models for reconfigurable intelligent surfaces: From surface electromagnetics to wireless networks optimization Proc. IEEE 110 2022 1164 1209 10.1109/JPROC.2022.3195536
42 Ulyanov D. Vedaldi A. Lempitsky V.S. Instance normalization: The missing ingredient for fast stylization Preprint at: arXiv 2016 10.48550/arXiv.1607.08022
43 Agarap A.F. Deep learning using rectified linear units (relu) Preprint at: arXiv 2018 10.48550/arXiv.1803.08375
44 Kennedy J. Eberhart R. Particle swarm optimization Proceedings of ICNN’95 - International Conference on Neural Networks 4 1995 1942 1948 10.1109/ICNN.1995.488968
45 Mirjalili S. Lewis A. S-shaped versus v-shaped transfer functions for binary particle swarm optimization Swarm Evol. Comput. 9 2013 1 14 10.1016/j.swevo.2012.09.002
46 Jiang W. Zhang Y. Wang R. Comparative study on several pso algorithms The 26th Chinese Control and Decision Conference 2014 CCDC 1117 1119 10.1109/CCDC.2014.6852332
47 Kingma D. Ba J. Adam: A method for stochastic optimization International Conference on Learning Representations (ICLR) (San Diega, CA, USA) 2015
