The terahertz (THz, approximately 0.1–10 THz) band has long represented a gap in research compared to the microwave and infrared bands, historically known as the “THz gap.” In recent years, THz technology has developed rapidly in radar and target identification^[1], thermal radiation^[2], wireless communications^[3], hazardous materials detection^[4], and non-destructive testing^[5]. THz absorbers are an important research area and are vital for applications such as target stealth and sensing technology^[6-8]. Due to the increasing complexity of electromagnetic (EM) interference in the real world today, it is particularly important that the absorption of the specified single-band can be done efficiently.

With the rapid development and advancement of micro/nanofabrication technology, metasurfaces such as artificial periodic subwavelength composite structures offer an ideal path to solve this problem^[9]. Numerous studies have shown that the combination of metasurfaces with different fine structures can achieve various absorption modes for THz waves, including single-band^[10], dual-band^[11], multi-band^[12], and broadband^[13-15]. However, in conventional design processes, researchers cannot anticipate the specific absorption indicators (e.g.,  absorption frequency, absorptance, bandwidth) of the designed absorber. They can only gradually find them in a continuous optimization process, which often requires tens or even hundreds of thousands of iterations of optimization to achieve. The process of optimizing the absorber parameters consumes time and money, and how to design for the required absorption frequency is one of the most pressing challenges to be solved.

The artificial neural networks (ANNs), as one of the most influential technologies of the last decade, can assist researchers in the design of a wide range of devices^[16,17]. ANNs typically consist of an input layer and an output layer, and one or more hidden layers in between them; the hidden layers can learn autonomously the relationship between the device structure and its optical counterparts (e.g.,  the complex impedance matching process in the absorption of EM waves and the excitation of electric and magnetic fields)^[18]. This enables researchers to rapidly determine the design of structures, simplifying the design process and increasing design efficiency^[19,20]. However, the current ANN-aided design of metasurface devices is limited to the design and optimization of the device^[21,22], with little research attempting to inversely design the device in conjunction with its requirements to achieve on-demand customization.

In this paper, we propose a conceptual method for the inverse design of the THz multifunctional metasurface absorber (MMA) based on the bidirectional deep neural network (BDNN). The BDNN (tandem neural network) consists of a forward modeling network and an inverse design network, which is able to predict the structural parameters of the metasurface with high accuracy, while also effectively solving the accuracy problem in the inverse design process^[23,24]. By employing the BDNN model, it becomes possible to achieve near-perfect absorption (absorptance greater than 99%) at the desired absorption frequency. As a proof of concept, one frequency of a relatively simple single-band metasurface absorber was selected for the inverse design, and the output parameters are validated by numerical simulations (midrange frequency, 5.0 THz). The same BDNN inverse design was used to superimpose the structure of the broadband metasurface absorber in the base of the single-band metasurface absorber to make the MMA. In addition, the absorption mechanism in different absorption modes is analyzed in conjunction with the electric field and effective medium theory (EMT), while the absorption spectrum of the MMA is analyzed at different angles of incidence. Unlike previous designs using the ANN to optimize parameters, this study enables the customization of the required MMA to suit the requirements, opening a new path for the design of high precision devices in combination with BDNN.

Figure 1(a) depicts the schematic diagram and side view of the structural unit of MMA, which adopts a simple three-layer structure design of graphene, polyethylene cyclic olefin copolymer (TOPAS), and silver (Ag) substrate, arranged in a top-to-bottom configuration. In the diagram, d1, d2, d3, and d4 indicate the length of the graphene square, the inner and outer edge lengths of the small graphene ring, and the outer edge length of the large graphene ring, respectively. The small and large graphene rings are divided by a thin strip with side length r, P indicates the period of the unit cell, and t1 and t2 denote the thicknesses of the TOPAS layer and the Ag substrate, respectively. All of the graphene except the small internal graphene squares is connected with thin graphene strips of 10 nm width to ensure continuity.

As a conceptual demonstration, the single-band absorption mode was first inversely designed using BDNN. Figure 1(b) presents the schematic diagram of an inverse design process based on a BDNN for generating the required parameters of the single-band metasurface absorber. The whole process consists of steps that have been described in detail in the Supplementary Material.

Graphene is a two-dimensional material, and its electrical conductivity can be determined using the classical Kubo model^[25]: $$\begin{gathered} σ(ω,μc,Γ,T)=je2(ω−j2Γ)πℏ2 \\ ·(1(ω−j2Γ)2∫0∞ε(∂fd(ε)∂ε−∂fd(−ε)∂ε)dε−∫0∞fd(−ε)−fd(ε)(ω−j2Γ)2−4(ε/ℏ)2dε), \end{gathered}$$(1)fd(ξ,μc,T)=(e(ξ−μc)/κBT+1)−1,(2)where ω, μc, and T represent the angular frequency, chemical potential, and temperature, respectively. Γ denotes the scattering rate, and fd (ξ, μc, T) represents the Fermi–Dirac distribution. The chemical potential (Fermi energy level) of graphene can be modulated by the application of a bias voltage. Figures 1(c) and 1(d) give the schematic of the voltage applied to the MMA and the calculated chemical potential as a function of voltage. The chemical potential of the internal and external graphene is controlled by two bias voltages, and each structural unit of the MMA is arranged in this way. The relationship between the chemical potential and the bias voltage is μc≈ℏvfπεrε0(Vg+V0)eh,(3)where V0 and Vg represent the natural doping-induced voltage offset (about −0.8V)^[26] and the applied bias voltage, respectively. Furthermore, the refractive index of the TOPAS material is 2.35 in the THz region^[27]. The relative permittivity of Ag in the THz region can be described using the Drude mode, as shown below^[28,29]: ε=1−ωp2ω(ω+iγ),(4)where ωp indicates the plasma frequency with the value of 1.37036×1016rad/s, and γ indicates the damping constant with a value of 1.71732×1016rad/s.

The fabrication process of the designed MMA was investigated in the following steps: (1) deposition of silver on a silicon substrate by magnetron sputtering to form the Ag substrate; (2) evaporation of TOPAS material by spin coating or electron beam evaporation on top of the Ag substrate to form the TOPAS layer; (3) growth of a graphene layer by chemical vapor deposition (CVD), which is then transferred on top of the TOPAS layer; (4) for patterning of the graphene, sputtering zinc onto the uppermost layer of graphene to remove the carbon layer selectively from those areas; (5) cleaning to obtain the proposed MMA, and then completing the full manufacturing process^[30].

In this study, the FNN was constructed using four units as inputs, representing the four parameters, while the two units corresponding to the obtained absorption spectrum were used as outputs [Fig. 1(e)]. The entire network comprised four fully connected hidden layers with 20, 200, 200, and 20 neurons, respectively. For the PNN, the two required absorption spectrum parameters were used as input units, and the corresponding parameters were predicted. The PNN also consisted of four fully connected hidden layers with 20, 200, 200, and 20 neurons, respectively. During the inverse design process, the PNN efficiently predicted the required parameters in approximately 30 seconds. This stands in stark contrast to the conventional design process, which typically requires weeks or even months for parameter selection. The utilization of the BDNN significantly accelerates the design process, making it more practical and effective.

During the training process, the learning rate (LR) of the BDNN was chosen as three values: 0.0001, 0.0002, and 0.0003, with each LR tested for 10,000 iterations. Figure S1 in the Supplementary Material displays the loss function and the coefficient of determination (R2) for both the training and test sets with different LR values. Remarkably, after 10,000 iterations, the BDNN achieved a significantly lower loss of approximately 10−3 when using an LR of 0.0002, outperforming the other two cases, which had a loss of around 10−2. In Fig. 2(a), a detailed comparison of the loss function and R2 between the training and test sets, employing an LR of 0.0002, is depicted. The results demonstrate that the BDNN exhibits a small generalization error and remarkable fitting capability, indicating its robust performance in the prediction tasks.

Figures 2(b) and 2(c) illustrate the outcomes of the predictions for the frequency and absorptance of the absorption peaks in the training set. In the analysis of the results, three common metrics were introduced to assess the predicted performance of the model: root mean square error (RMSE), mean absolute percentage error (MAPE), and standard deviation (STD), respectively (see the Supplementary Material). The RMSEs for the predicted absorption peak frequency and absorptance in the test set were only 0.0542 and 0.0051, with MAPEs of 0.28% and 0.40%, respectively. The quantitative analysis conducted demonstrates that the BDNN exhibits remarkably high prediction accuracy for MMA prediction and this prediction method is feasible and effective.

Subsequently, the inverse design was conducted using the BDNN to determine the frequency and absorptance of various absorption peaks. To simulate the realistic needs of the working environment, the inverse design was performed for 10 frequency points within the range of 0.5–10 THz, with intervals of 1 THz (from 1 THz to 10 THz). While predicting the absorption peak frequencies, the absorption peaks were required to have an absorptance of more than 0.99. The parameters corresponding to each absorption frequency, as predicted from the BDNN inverse design, are detailed in Table S2 in the Supplementary Material, and the absorption curves at different frequencies obtained from the simulations based on the predicted parameters are given in Fig. 2(d). Figures 2(e) and 2(f) show the errors between the predicted absorption spectra, concerning absorption peak frequencies and absorptance, respectively, and their expected values. Notably, the absorption peak frequencies acquired from the BDNN inverse design demonstrate a high degree of correspondence with the targeted frequencies, with very minimal errors, and the RMSE and MAPE values are 0.0037 and 0.09%, respectively. Moreover, the single-band metasurface absorber resulting from the BDNN inverse design achieves an absorptance of 0.9978, with an STD of 0.0016, exhibiting excellent stability while ensuring exceedingly high absorptance. The single-band metasurface absorber obtained through the BDNN inverse design covers a broad range of frequencies, reaching up to 9.5 THz, covering the quasi-entire THz gap.

To demonstrate the conceptual advantages of using the BDNN for the inverse design of parameters, it is necessary to compare it with some other classical algorithmic models. Two classical algorithms are selected for testing, the K-nearest neighbor (KNN) and the random forest (RF) algorithms, which have been extensively utilized in various applications, including spectral identification and classification, solar spectrum prediction, and water resource management ^[31,32]. As with the BDNN, the KNN and RF algorithm models were inversely designed using the same original data set, and the training set was scaled to 0.9. In addition, both algorithm models were hyper-parametrically tuned using a grid search to ensure optimal performance prior to the inverse design.

Similarly, the frequency and absorptance of the absorption peaks in the test set were tested using KNN and RF algorithms prior to the inverse design, as shown in Figs. S2(a)–S2(d) in the Supplementary Material. The RMSEs predicted by the two algorithms for the absorption peak frequencies were 0.3514 and 0.0751, and the MAPEs were 4.13% and 1.47%, respectively. For the prediction of absorptance, the RMSEs for the two algorithms were 0.0594 and 0.0434, with MAPEs of 4.31% and 3.15%, respectively. In comparison to BDNN, the KNN and RF models exhibited larger errors and lower prediction accuracy, which also proved the conceptual advantage of the BDNN. In addition, Fig. S2(e) in the Supplementary Material presents the error statistics for absorption peak frequency and absorptance for all test sets (100 times for each parameter), and it is evident that the BDNN predicted results have a low error, with the number of errors less than 1% reaching 83 and 73, respectively, achieving remarkably high accuracy. And 50 repetitions of each algorithm were performed to verify stability, with the BDNN demonstrating extremely high stability while maintaining a high R2 [Fig. S2(f) in the Supplementary Material].

Subsequently, the KNN and RF models were employed for the inverse design of absorption peaks at different frequencies, and the parameters obtained are presented in Tables S3 and S4 in the Supplementary Material, respectively. Figure 3 showcases the absorption spectra obtained from simulations based on the parameters derived from the inverse design, with both algorithm models showing varying degrees of error. Further quantification of errors in frequency and absorptance of the absorption peaks obtained from the inverse design of the KNN and RF models is illustrated in Fig. 3, yielding RMSE values of 0.1691 and 0.1527 for frequency errors, MAPE values of 4.68% and 3.48%, and STDs of 0.0118 and 0.0049 for absorptance errors, respectively. Compared to KNN and RF, the BDNN is able to maintain high accuracy in simultaneous two-parameter prediction, demonstrating a great advantage.

After theoretically confirming the feasibility of the BDNN inverse design, the method is applied to the design of the broadband metasurface absorber. Based on the structural parameters at 5.0 THz and keeping all the structural parameters unchanged, the structural design of the broadband metasurface absorber is carried out using the external space of square graphene. A simple double square graphene ring is selected for the structural design. Similar to the previous inverse design of the structural parameters using the BDNN, in the inverse design of the structural parameters, their value ranges and minimum accuracies are shown in Table S5 in the Supplementary Material, respectively. The difference from the single-band metasurface absorber is that, in addition to the absorption peak frequency and absorptance, the effective absorption bandwidth (absorptance greater than 90%) is added as the target, and it is expected to obtain the broadband metasurface absorber with 5.0 THz as the center frequency and sufficiently large absorption bandwidth. It is worth noting that the chemical potential of graphene is no longer taken as a variable parameter; instead, the external square ring graphene is uniformly set to 1.0 eV, and the internal square graphene is set to 0 eV. In other words, the MMA is in the single-band absorption mode when the chemical potential of internal and external graphene is set to 0.39 eV and 0 eV, respectively, and in the broadband absorption mode when the chemical potential of internal and external graphene is set to 0 eV and 1.0 eV, respectively. Similarly, the optimal parameters obtained from the inverse design are d2=3.97μm, d3=7.95μm, and d4=14.76μm. Numerical calculations using the output structural parameters verify that the broadband metasurface absorber can achieve effective absorption in the range of 3.18–6.96 THz with an absorption bandwidth of 3.78 THz.

To verify that the structural parameters are at their optimal values, each parameter is verified separately. Figures S3 and S4 in the Supplementary Material show the absorption spectra calculated in single-band and broadband absorption modes, respectively, for each different structural parameter. A quantitative comparison of the statistically obtained optical properties at different parameters is given in Fig. 4; in single-band absorption mode, the absorption rate and frequency are calculated separately for different parameters, and in broadband absorption mode, the absorption bandwidth is additionally calculated. It should be noted that, in the validation of each parameter individually, the other parameters are taken as predicted and kept constant. In single-band absorption mode, the absorption is at its maximum value (99.61%) for all structural parameters, and the absorption frequency is located at 5.0 THz. In the broadband absorption mode, the effect of d2 on the single-band absorption mode is larger because the small distance affects the square graphene surface plasmon resonance (GSPR), while the changes of d3 and d4 do not affect single-band absorption mode. All the structural parameters are also in the optimal values, which verifies the great advantage of the BDNN in inversely designing optical micro-nano devices.

To gain insight into the absorption process of THz waves by the MMA, an analysis of the absorption mechanism is conducted in conjunction with the EMT. In a complex system, the transport of waves can be approximated and described by the EMT. The normalized equivalent surface impedance Z of the MMA can be expressed as follows^[33]: Z=±(1+S11)2−S212(1−S11)2−S212,(5)R(ω)=((Z−Z0)/(Z+Z0))2,(6)where S11 and S21 represent the reflection and transmission coefficients of the incident THz wave and Z0 denotes the impedance of free space. For impedance analysis, the absorption spectra of the absorption frequency point (5.0 THz) were selected, and Figs. 5(a)–5(c) present the calculated absorption spectrum and normalized equivalent surface impedance in both single-band and broadband absorption modes. It is evident that at the absorption frequency, the real part (green solid line) and the imaginary part (yellow dashed line) of the normalized equivalent surface impedance converge to 1 and 0, respectively, indicating the impedance match state between the MMA and free space. Consequently, the incident THz wave can enter the MMA without any reflection, while the metal substrate prevents the transmission of the THz wave, achieving perfect absorption of the incident wave.

Besides, to investigate the absorption mechanism of MMA in more depth, the calculated normalized electric field at the absorption peaks of the selected frequency absorption spectra was plotted, as shown in Figs. 5(d)–5(i). The electric field distribution in the x–y plane shows that the excited electric field in single-band absorption mode is mainly located at the edges of the square graphene, with particularly pronounced electric fields observed at the four top corners of the graphene, which is a typical feature of the GSPR. In broadband absorption mode, the excited electric field is mainly located at the gap of the external graphene ring, which indicates that the broadband absorption of MMA originates from the GSPR of the external graphene. Additionally, the x–z cross-sectional analysis of the electric field reveals that the energy loss primarily occurs at the surface of the MMA graphene layer^[34,35].

The impact of different azimuthal angles on the vertical incidence of THz waves is also important, as depicted in Figs. S5(a) and S5(d) in the Supplementary Material. Due to the perfect symmetry of the designed structure, the MMA maintains a stable operating state at any angle of polarization incidence without being affected, exhibiting remarkable polarization insensitivity. In addition, the effect of different incident angles on the absorption performance of the MMA in TE- and TM-polarization modes has been investigated [see Figs. S5(b), S5(c) and Figs. S5(e), S5(f) in the Supplementary Material]. At small incidence angles (less than 40°), the MMA consistently maintains an absorption above 0.9 with stable absorption performance. As the angle gradually increases, the absorptance of MMA decreases, and the absorption spectra show different degrees of blue-shift in both single-band and broadband absorption modes (especially in the TM-polarization mode). This is because, at large incident angles, parasitic resonances are generated in the graphene layer, and the resonant frequencies are shifted toward high frequencies^[36]. In conclusion, the MMA obtained by the inverse design can meet the needs of wide-angle incidence in practical applications.

Finally, we try to combine BDNN for dual-frequency absorber customization. In performing the inverse design of the two absorption peaks, a second square graphene was assembled at the edge of each structural unit in addition to the square graphene at the center of the metasurface, as shown in Fig. S6(a) in the Supplementary Material. Taking 4.0 THz and 5.0 THz as examples for the inverse design, the obtained absorption spectra are shown in Fig. S6(b) in the Supplementary Material, and the center frequencies of the absorption peaks are located at 3.989 THz and 4.998 THz, respectively, with a very high prediction accuracy. This also proves that it is feasible and effective to combine the BDNN with the metasurface for the inverse design of multi-frequency absorbers.

In conclusion, we propose a novel method based on the BDNN that facilitates the accurate and on-demand inverse design of MMA to achieve multifunctional absorption at specific frequencies. The BDNN, obtained by concatenating the FNN and PNN, can accurately predict the parameters of the MMA. Ten frequency points in the range 0.5–10 THz were verified by simulation, which showed that the RMSE and MAPE of the frequency errors were only 0.0037 and 0.09% and that the average absorptance reached 0.9978. On this basis, a frequency point was selected for the inverse design of MMA using the BDNN, and the absorption mechanism was analyzed. The MMA obtained by the inverse design can achieve single-frequency perfect absorption and at the same time can switch to broadband absorption with an effective absorption bandwidth of 3.78 THz. At different angles of incidence, the MMA exhibited polarization insensitivity and wide-angle incidence. These findings demonstrate the feasibility and effectiveness of employing the BDNN for the on-demand inverse design of metasurfaces, offering a promising avenue for the design of ultra-compact precision multifunctional instruments with exceptionally high accuracy.