Frequency transfer and inverse design for metasurface under multi-physics coupling by Euler latent dynamic and data-analytical regularizations

Domain: earth_science
Published: 2025-03-06
Editorial outcome: major_revision

Abstract
Frequency transfer is a key challenge in machine learning as it allows researchers to go beyond in-range analyses of spectrum properties towards out-of-the-range predictions. Traditionally, to predict properties at a specific frequency, targeted spectrum is included in training data for a deep neural network (DNN). However, due to limitations of measurement or computation source, training data at some frequencies are hardly accessible, especially for multi-physics problems. In this work, we propose a multi-physics deep learning framework (MDLF) consisting of a multi-fidelity DeepONet, a Euler latent dynamic network, and a data-analytical inversion network. Without the knowledge about multi-physics response, MDLF is successfully generalized to unseen frequency bands for both parametric and free-form metasurface by dynamically utilizing a Euler latent space and single-physics information. Moreover, an inversion method is introduced to incorporate hybrid a priori in inverse design of metasurface. Under EM-thermal coupling, we verify the proposed MDLF numerically and experimentally.

Introduction
Nowadays, artificial intelligence (AI) have attracted many interests in bioscience1, material science2, and energy science3, which improves the quality of detection and analysis4, accelerates new discoveries5,6, and achieves real-time control7. In particular, as a surrogate to full-wave simulation8–10, the data-driven artificial neural networks11–13 (ANNs) have shown extraordinary potential in speeding up design procedure14,15 and shortening the gap between measurement and design16,17, which greatly relieves the pressure of numerical simulations to obtain a suitable solution in inverse design18–22.

For example, variational auto-encoders (VAE)23,24 is thought to be a practical pre-processing method to compress the data space for better searching the design space by another optimization method. Developed from the generative adversarial network (GAN)25, conditional GAN (cGAN)26, as one of the most popular frames in inverse design these years, constructs both forward and inverse surrogate models by adversarial strategy, partly alleviating the problem of non-uniqueness.

In electromagnetic (EM) society27–29, a traditional ANN-based surrogate model builds a mapping from the material distribution of objects to the spectrum properties, in which the prediction range is determined by the training data and is fixed once the network is trained. However, due to the limitations of measurement or computation source, the training data sometimes are not complete to cover the scope of interest, and the performance of deep learning deteriorates dramatically when it comes to out-of-the-range tests, especially for complex systems or models, such as multi-physics coupling problems. Therefore, generalizations of the surrogate model, which allow researchers to go beyond in-range analyses of properties toward out-of-the-range predictions, are crucial.

Due to the unique ability to manipulate EM waves, metasurface has shown its promising abilities in hologram30,31, absorber32, semiconductors33, 5G/6G base station34–36, cloaking37, radome38 etc. Generalizations with respect to the operating state of the metasurface, like polarization transfer and dispersion transfer39, have attracted many attentions to relieve the burden of accumulating dataset and improve the applicability of the model. Similarly, frequency transfer, which is to predict the unknown spectrum given part of the frequency data, could be divided into two types according to the frequency range provided by the training dataset.

a. In the first type, the range of training datasets already covers the predicted frequency: For example, the mapping relationship between the low-frequency and high-frequency wave field is built in Geophysics fields40, where high-frequency data as the output is included in training data. Furthermore, in the field of metasurface, a data-driven spectra-to-spectra strategy is introduced to decode the unknown spectrum of metasurface from the latent space with the known spectrum41. By operator learning mechanism, the interpolation of frequency is achieved42.

b. The second type of frequency transfer is out-of-distribution (OOD) generalization, i.e., the predicted spectrum is not seen in the training stage, which is the problem discussed in this work. There are two mainstream approaches, i.e., multi-fidelity modeling and fine-tuning43. For example, transfer learning is adopted in fine-tuning the backbone to achieve OOD with a sparse amount of data given in the extrapolation44. In addition, there are also physics-based approaches to extrapolation. Based on the correlation of metasurface size scaling and resonant wavelength, the generalizability in operation wavelength in the forward method is achieved39. Based on a priori provided by the equivalent circuit, spectral extrapolation of metamaterials is achieved45.

Despite the convincing achievement in the above literature, the mentioned works only focus on a single physical response. However, with the decreasing sizes and increasing power of electronic devices, forward modeling and inverse design of metasurface under multi-physics coupling46–48 becomes demanding. In the system where the EM-thermal effects dominate47, for example, as the power consumption increases, the temperature rise effects are unavoidable, which leads to the change of material constitutive parameters and further makes EM response deviate from the single-EM simulation, i.e., numerical calculation of object/system’s EM behaviors governed by Maxwell’s equations. However, modeling and designing metasurface considering both the EM and thermal multi-physical properties, i.e., multi-physical simulation, are challenging due to the complicated physical process and high computational costs. Specifically, to include EM propagation, heat loss, and heat transfer in the metasurface, Maxwell’s equations and heat conduction equation are alternately solved in the steady-state solution of the EM-thermal coupling49, between which temperature-varying parameters build the bridge.

Considering the difficulty in accumulating high-quality datasets, the key to a multi-physical surrogate model with high precision as well as frequency transfer is improving the data efficiency in the context of scientific machine learning (SciML). On the one hand, by explicitly introducing the governing equations into the loss function, physics-informed neural networks (PINNs)50–52 have made good practices in solving PDEs given specific questions. But, in this case, the multi-physical behaviors in the metasurface are governed by different partial differential equations (PDEs), which is hard to decouple. On the other hand, by learning operator mappings between functions, neural operators, such as deep operator networks (DeepONet)53–55, and Fourier neural operators (FNO)56, have shown tremendous success in modeling complex physical phenomena and engineering problem11,57,58, where the capability to achieve parameter efficiency59,60, small data set training61–63, and to capture non-local interactions64 are witnessed. Compared with traditional ANNs, it is reported that DeepONet could learn operators accurately and efficiently from a relatively small dataset in practice61. Also, to alleviate the pressure of high-quality data accumulation, a multi-fidelity strategy is adopted62,63. In predicting the characteristics of time-varying systems, latent dynamics networks are developed to discover the system evolution of spatio-temporal processes, where non-Markovian effects are captured by tracking the system history64.

In this work, we propose a multi-physics deep learning framework (MDLF), as shown in Fig. 1. The forward networks consist of two parts, i.e., multi-fidelity DeepONet to predict EM S-parameter performance (S M) and a latent-dynamic network to predict temperature performance (T(f)), where both have considered EM-thermal multi-physics coupling effects. More specifically, since the computational cost is high when accumulating training data for predicting S M, the multi-fidelity DeepONet65 is used to reduce the need of high-fidelity data (S M), where data of S-parameter performance without EM-thermal coupling effects (S E) are firstly accumulated to train a low-fidelity network, i.e., Net 1 in Fig. 1. The output of the low-fidelity network is further input into the high-fidelity network (Net 2) and formed as a residual part to reduce the difficulty in training Net 2. Further, since obtaining thermal effects at high frequency is either computational or experimental challenging, we further developed a latent-dynamic network64 (Net 3 & Net 4) to predict temperature performance (T(f)) considering EM-thermal coupling effects. By using an Euler-method-based dynamic latent, the proposed network can predict temperature performance at a higher frequency not seen in training data, i.e., out-of-distribution (OOD) prediction. In the inverse method, since designing a metasurface considering both EM and thermal functionalities is a challenging multi-objective problem, hybrid regularizations66 are introduced in the inverse network (Net 5) to boost the training process. More specifically, besides the loss defined on geometry parameters, forward solvers, including EM and thermal predictor and an analytical model, are further introduced to alleviate non-unique problems and constrain resonant points, respectively.

Results
Network architecture
We begin with a brief introduction of the vanilla DeepONet61 and multi-fidelity DeepONet65 proposed by Lu. et al. and describe some architectural choices in this work. Different from the traditional ANNs, a DeepONet consists of two sub-networks, i.e., branch net and trunk net, where the branch net encodes the input function at a fixed number and trunk net encodes the locations for the output functions61. Given the input function v and the locations ξ, the output function can be approximated by merging the output of branch net and trunk net:
G(v)(ξ) ≈ ∑ b k t k (1)
k

At the same time, to make full use of both low-fidelity and high-fidelity data, the multi-fidelity DeepONet is proposed with low-fidelity and high-fidelity models, in which it is recommended to append evaluation of low-fidelity in the trunk-net inputs and leave the branch-net unchanged, i.e., ref. 65,
G H (v)(ξ) = Q(v)(ξ, G L (v)(ξ)) (2)
where G H represents the high-fidelity model, G L represents the low-fidelity model.

In this case, by taking frequency (f) as ξ in the trunk net and geometric input (G) as v in the branch net, DeepONet is a suitable frame to achieve a free sampling of frequency and a continuous representation of the spectrum in the metasurface. As shown in Fig. 2a, to predict S M, the multi-fidelity DeepONet contains one low-fidelity network, Net 1, and one high-fidelity network Net 2, which are trained separately. Firstly, with the framework of unstacked DeepONet, the information from the geometry domain G of the multi-physic metasurface and frequency domain f are processed through two independent channels in Net 1, and the S-parameter by single-EM simulation S E is predicted.

Then, the output of the trained Net 1 is further input into the high-fidelity network (Net 2) and formed as a residual part (S R) to reduce the difficulty in training Net 2. In the trunk net, the EM response by Net 1 corresponding to the current f is added, as well as that by two frequency points before and after the current frequency indicating the local features. Considering low fidelity Net 1 provides a preliminary knowledge of space mapping and Net 2 only learns a small value, it allows a sparse sampling in the training set of Net 2, which reduces the pressure of multi-physical simulation. Moreover, in this work, the same geometric input (G) goes over the progressive feature extraction since the feature of G is not fully extracted in Net 1, and there are totally different S R with the same input of S E and f, but different G, as shown in Supplementary Fig. 5.

In predicting the temperature performance (T(f)) under EM-thermal coupling, latent dynamics networks64 are developed in the frequency sequence flow. As presented in Fig. 2c, a(f) is computed from the parameters in the latent space, where no label of a(f) is provided in the whole training process, and the network automatically discovers a compact representation of the learning task, i.e., latent space. More specifically, similar with the frame of the encoder and decoder, Net 3 in Fig. 2c encodes S E into latent space, and Net 4 decodes a(f) into the temperature output T(f). In this way, the weights of Net 3 are updated through the back propagation of loss between the outputs of Net 4 and its true value. Besides, given the initial state, i.e., a(f) = 0, the update of latent space in frequency sequence follow the Euler method,
a(f n+1) = a(f n) + ∂a/∂f × β × Δf, Δf = f n+1 - f n (3)
where Δf represents the step in the frequency sequence, β is the hyperparameter for regularization, and ∂a/∂f is the changing rate of latent space, which is also the output to Net 3.

To conclude, given no label for the latent space, Net 3 specifies the update rules of a(f) in the frequency sequence, Net 4 ensures the effectiveness of a(f) in predicting the temperature at each frequency point, which together achieve frequency transfer to the unknown frequency. Besides, in Net 4, we also use S E as one of the inputs because the temperature results and the EM results are highly correlated. However, S E is not the only determining factor to T(f). For instance, even with the same EM response at different frequencies, structure can appear in different thermal behaviors, as shown in Fig. 3c.

Forward Method
Firstly, we consider parameterized metasurface structures, whose predefined topology is given in Supplementary Fig. 1. In total, with 720 individual geometries, low-fidelity data S E are obtained with a sample interval of 0.1 GHz within 2–14 GHz, high-fidelity data S M are obtained with a sample interval of 0.2 GHz, which are randomly divided into training set and validation set. As shown in Fig. 3a, for a typical geometry in the trainingset, a comparison between the S M by multi-physical full wave simulation and those by multi-fidelity model and low-fidelity model is given, which indicates the differences between S E and S M. For the test out of the training set, the comparison between the full-wave simulation and the predicted S M is given in Fig. 3b, which matches well. It is noted that the sample interval in testing the multi-fidelity model is 0.1 GHz, which means frequency interpolation in S M.

To quantify the overall accuracy for predicting EM response, relative error is used, i.e.,
E R = ∑ |y i - f(x i)| / ∑ |y i| (5)
where y i represents the true value and f(x i) is the predicted value. Tested with 106 random samples unseen in the training process, the multi-fidelity DeepONet obtains an average relative error of 0.036. Then, the network trained only with high-fidelity data (single-fidelity DeepONet) is compared with the proposed multi-fidelity model to illustrate the importance of the proposed multi-fidelity DeepONet. More specifically, the trained single-fidelity DeepONet is tested with the same validation set with 106 random samples, where the relative error is 0.154.

By adjusting the sampling strategy of high-fidelity data, the influence of sampling interval ratio between low-fidelity model Net 1 and high-fidelity model Net 2 is studied. As indicated in the left of Table 1, as the sampling interval ratio increases from 2:1 to 12:1, the relative error in the validation set rises like a step. For example, the error in Case-R2 and Case-R3 is very close, however, from Case-R1 to Case-R2 and to Case-R4, the rise in error is more obvious, since their interval ratio increases exponentially. A similar situation can be found in Case-R3 and Case-R6. In addition, if the training data is randomly selected in the original high-fidelity data, which has 61 points in one single geometry, the influences of the sampling number in high-fidelity data are shown on the right of Table 1. Therefore, with a preliminary knowledge of space mapping provided by the low-fidelity model Net 1, the multi-fidelity strategy allows a relatively smaller high-fidelity dataset.

Average temperature by latent dynamics networks
To illustrate the frequency transfer achieved by latent dynamics networks, we remove the data whose corresponding frequency is over 10 GHz from the training set, that is, the temperature prediction from 10.2 to 14 GHz is from frequency transfer. Similarly, relative error in Eq. (5) is used to quantify the overall accuracy for predicting average temperature. In the trained networks, E R for the temperature within 2 to 10 GHz is about 0.007, and for the temperature in 10.2 to 14 GHz is about 0.043.

As shown on the left of Fig. 3c, the comparisons between the average temperature of the geometry by full-wave simulation and prediction are given, where EM response as the inputs to Net 3 and Net 4 are also presented. From the figures, the extreme points of temperature have a high correlation with the extreme points of EM response in the frequency domain. Yet, even with the same EM response, geometry appears at different temperatures at different frequencies. For frequency transfer, in figures, the areas with background color are totally unseen in training, in which the simulation and prediction matches well. Besides, it is noted that the most dramatic temperature effect occurs near the extreme point.

To evaluate the contribution of latent space a(f) in temperature prediction as well as frequency transfer, a contrast test between the traditional method without the proposed latent space and the proposed method with latent space is conducted. In the traditional method, f, EM response, and G are the inputs to DNN and T(f) is the output, which has a similar structure with Net 2 but has a different output. Specifically, in the traditional method, with the same dataset, E R within 2–10 GHz is about 0.0281, and in 10.2–14 GHz is about 0.1643, which is much larger than the error in the proposed method.

To further quantify the prediction error at each frequency point, square error is used. As shown in Fig. 3d, in the proposed model, part of the structure samples is used for training, which means out-of-range tests are carried out in both geometry and frequency range. Compared with the results from a traditional DNN method without latent space, the proposed method obtains a smaller test error. In addition, a typical example of temperature prediction is given in Fig. 3e, where the comparisons between the full-wave simulation, the proposed method, and the traditional method are conducted. From the figure, even if the curve by the traditional method meets the trend of the simulated result, the predicted temperature deviates from the simulation, when it comes to the extrapolated frequency points.

Furthermore, besides the parameterized metasurface structures, by changing full connected network to a convolutional neural network to deal with 2D patterns while keeping the framework of the prosed method unchanged, the frequency transfer of the proposed method on free-form meta-atoms with a higher degree of freedom are also tested. For example, the thermal performance of free-form single-layer meta-atoms are predicted, as shown in the middle of Fig. 3c, where the extrapolated frequency band is set to be 11.2–14 GHz with a relative error of 0.047. For more complex EM-thermal coupling effects, a three-layer meta-atom is studied in this work. In this case, the top and bottom layers are the same 2D pattern with a high degree of freedom, as shown on the right of Fig. 3c, where the simulation and prediction match well. The results from other geometries are given in Supplementary Fig. 8. From these experiments, it is shown that the frequency transfer of the proposed method still validates in higher dimension cases.

Inverse design by data-analytical networks
In the development of an inverse design network, the weights in Eq. (4) are important to direct the training. Under the same training epochs, as shown in Fig. 4a and b, the relative errors in the validation set are given under different training cases, in which the numbers from left to right in parentheses represent the weights of data-driven regularization w 1, analytical model regularization w 2, EM forward model regularization w 3, and temperature forward model regularizations w 4, respectively. Besides, E1 represents the relative error on geometry determined by data-driven regularization. E2 represents the relative error on the resonant frequencies determined by analytical-driven regularization, in which E21 and E22 represent the relative error on the first and second resonant frequencies, respectively. E3 represents the relative error on EM response predicted by the EM forward model. E4 represents the relative error on average temperature predicted by the thermal forward model.

As shown in Fig. 4b, compared with Case 3, driven by data and EM forward model, Case 2, which is driven by data and analytical model, has a similar E21 and a lower E22, which means a better prediction in resonant frequency. From Case 1 to Case 3 and Case 5, the EM forward model and analytical model are added into the traditional method in turn, in which a gradual accuracy improvement can be observed. In addition, Table 2 gives the detailed value in Fig. 4. From the table, Case 5, which is adopted as the proposed data-analytical-driven model, performs the best in many indexes and has a very close value with the best performance in other index.

Furthermore, the effect of analytical-driven regularization is visualized in Fig. 4c and d. As shown in Fig. 4c, although the EM response in the passband obtained by Case 2 (governed by the analytical model) is less accurate than that by Case 3 (governed by the EM forward model), a higher accuracy in pole prediction is witnessed in Case 2. From Fig. 4d, by adding analytical model to the model driven by data and EM forward method, the inverse design achieves more accurate local performance deign while ensuring the overall consistency.

In the inverse design scenario, usually, we do not have detailed performance distributions, i.e., desired EM performance, before we design. Instead, we usually have some incomplete requirement, for example, we want to find a geometry with two passbands at about 4 GHz and 8–10 GHz, thus a polyline of performance is draw according to these obscure requirements, as shown in Fig. 4e, indicating the desired passband location and bandwidth. Besides, the extra peak between the two passbands of the parameterized metasurface in Fig. 3c is eliminated for the purpose of achieving stability in EM and thermal performance. However, no one knows whether a metasurface structure satisfying the plotted polyline exists or not. In this case, the aim of inverse method here is to find a geometry to satisfy the fictitious performance as much as possible, deviations outside the required passband could be acceptable. There by, by applying Net 5, a comparison between the desired EM response and the simulated response from the designed geometry is given in Fig. 4e, where a small increase can be found at the end of frequency, considering the EM properties of the geometry itself.

Furthermore, in the situation where real EM response curves from the simulation are available, consistency between the desired EM response and the EM response from the designed geometry is observed, as shown in Fig. 4f. Particularly, the characteristic of a peak at the end of the frequency is captured by the inverse method, and the inversely designed geometry can reproduce this feature.

To further consider both EM and thermal multi-physics functionality in optimizing metasurface structures, the proposed model with Input-T, i.e., the value and position of maximum temperature within 2–8 GHz, as one of the inputs, is compared with the model trained without Input-T. In this case, by adding Input-T, we are targeted to constrain the thermal performance in the metasurface besides EM performance under EM-thermal coupling effects. As shown in Supplementary Fig. 10, under the same training weights, the relative errors in the model trained without Input-T are larger than those in the proposed model, especially in thermal performance. Specifically, with the same desired EM response, the achieved performances by the results from the inverse model trained without Input-T (Option 1) and with Input-T (Option 2) are given in Fig. 5a, b, respectively, whose design parameters are given in Supplementary Note 4. As shown in Fig. 5a, Options 1 and 2 have similar passband frequency characteristics, which are very close to the desired one. However, as shown in Fig. 5b, by constraining the position of the maximum temperature near 4 GHz, Option 2 performs better in the average temperature compared to Option 1, where the maximum temperature rise occurs at a higher frequency in Option 1 meaning more serious thermal effects happen.

Experimental results
To verify the thermal performance by the inverse method experimentally, design results in Fig. 5a, b are fabricated with 4*2 units, whose total size are 40 mm * 20 mm, as shown in Supplementary Fig. 13c. Then, the fabricated samples are irradiated in the waveguide with power amplifier (PA) enhancing the microwave from vector network analyzer (VNA), whose measurement platform is given in Supplementary Fig. 12c in Supplementary Note 5. After 3 mins of continuous wave irradiation at 7.9–8 GHz and 20 W, the surface temperature distributions of samples are given in Fig. 5c, where the highest and average temperature of Option 1 (347.45 K, 321.22 K) are higher than those of Option 2 (331.85 K, 311.71 K), while the environment temperature is close. In addition, the measured temperature at the quarter and the half of the samples (Line1, Line2) are compared as shown in Fig. 5d, where the temperature of Option 1 (design without Input-T) is generally above that of Option 2 (design with Input-T), with similar measured EM response as shown in Supplementary Fig. 14c.

To verify the accuracy of the proposed method, a parameterized geometry inversely designed by Net 5 is fabricated with 20*20 units, as shown in Supplementary Fig. 9 and 13a. Transmission (S 21) under low-power input (S E) is measured and compared with simulated S E, as shown in Fig. 6a.

In addition, the temperature rises of the fabricated sample under higher power EM radiation are measured, in which the input to the horn antenna is about 40 dBm, as shown in Supplementary Fig. 12b. For each single investigated frequency, the radiation time is about 20 to 30 mins when the temperature rises from room temperature to a steady state. The measured results are calculated by averaging the temperature rise in the whole surface of the sample, as shown in Fig. 6b, whose trend is close to that predicted by Net 3 and Net 4. The main reason for the difference in the numerical scale is that the input power to one single metasurface unit cannot reach the level in multi-physical simulation. At 3.65 GHz, as shown in Fig. 6c, when the thermal effect is most salient, the maximum temperature is in the center of the sample surface, indicating the maximum power received is at the same position. At 3.2 and 5.2 GHz, however, when the temperature rise is less apparent, the temperature distribution is relatively uniform, as shown in Supplementary Fig. 15.

The measurement of the EM response under high-power irradiation is conducted in the waveguide with PA. The inversely designed result is fabricated with 4*2 units with the size of 40 mm * 20 mm, as shown in Supplementary Fig. 13b. After 9 mins continuous wave irradiation at 6–7.05 GHz and 20 W, the measured S 21 are basically unchanged during irradiation, as shown in Fig. 6d since the heat effect in the design sample is designed to be not severe in the observed frequency, whose highest and average temperature after 10 mins irradiation are 314.15 K, 305.12 K. However, due to the difference in the EM behaviors of small samples in waveguides and periodic large arrays in free space, the measured S 21 in waveguide is not compared with that in Fig. 6a, which is the measurement in free space.

To further highlight the importance of Net 2, which is the importance of predicting the changes of EM response under multi-physical coupling, a contrast sample is further fabricated with the same size of the designed sample and measured in the same condition, as shown in Supplementary Fig. 13b. As shown on the left of Fig. 6e, after 10 mins irradiation under 20 W microwave of 6–7.05 GHz in waveguide, the surface temperature distribution of contrast sample is given, where the highest and average temperature is 360.15 K and 342.17 K. The increase of temperature further leads to the change of the EM properties of the sample. As shown on the right of Fig. 6e, it is obvious that S 21 of the contrast sample is shifted to the left as the irradiation time increases and basically stable at about 9 min, where S 21 at 0 min is S E and that at 9 mins is S M. Therefore, it is seen that EM response is significantly affected by multi-physical coupling and Net 2, that considers multi-physical coupling is important.

Methods
Metasurface under EM-thermal coupling
Apart from energy flow considered in the traditional works, where some of the input energy is reflected or passes through the structure into the free space, in this work, the energy absorbed by the lossy materials and conductors is also analyzed, which is transferred into heat, as discussed in Supplementary Note 1. To quantify EM-thermal coupling, the temperature-related material model is established. For example, the relative dielectric constant of the substrate is descried by the Debye model, as given in Supplementary Equation (1). The relationships between temperature and the conductivities (including both electrical conductivity σ Cu and thermal conductivity κ Cu) are fitted in the form of polynomials according to the experimental data67, as given in Supplementary Equation (2,3).

The unit structure of the parameterized metasurface contains three layers based on the substrate with a height of 4 mm, as shown in Supplementary Fig. 1. The metal patterns on the top and bottom layers are the same, and the pattern on the middle layer is a circular ring, where some of the structure parameters are fixed, and the other parameters are treated as design parameters, i.e., geometry (G).

The single-layer free-form meta-atoms in Fig. 3c are also based on the mentioned substrate with a height of 1 mm. The 2D random patterns are generated by Gaussian random fields15, where a center-symmetric pattern is produced by mirroring a one-eighth square pattern, as shown in Supplementary Fig. 3. The three-layer free-form meta-atoms in Fig. 3c are based on the mentioned substrate with the height of 4 mm, in which the top and bottom layers are the same 2D pattern and the middle layer is a circular ring, as shown in Supplementary Fig. 4.

Numerical simulation and training details
Both single-EM and EM-thermal simulations are carried out in the commercial software package COMSOL Multiphysics. For parameterized metasurface, design parameters are swept in the data collection, where 720 independent geometries are obtained. For single-EM data, frequency sampling from 2 to 14 GHz with an interval of 0.1 GHz (121 points) is adopted. For EM-thermal data, the sampling interval is 0.2 GHz (61 points), which costs about 7200 mins in total. On the contrary, the proposed learning-based solvers predict the results within 1 second. For free-form meta-atoms, COMSOL Multiphysics with MATLAB is used, where 544 independent single-layer meta-atoms and 648 independent three-layer meta-atoms are obtained. For the data pre-processing, the EM response is translated into the magnitude, and the average temperature is mapped by log10.

The multi-fidelity DeepONet Net 1 and Net 2 are developed with the DeepXDE library under Tensorflow 2.13. The latent dynamics networks Net 3 and Net 4 and the inverse design network Net 5 are developed under Keras and Tensorflow 2.13. The details of the architectures, as well as the network training, can be found in Supplementary Note 2, 3, and Supplementary Table 1. The time of training the networks are about 62 mins, where the loss decreases as the training process is given in Supplementary Fig. 6. It is worthwhile to mention that both accumulating dataset and training the networks can be done offline, and accumulating dataset can be further done in parallel.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The paper is not well written and hence, it is difficult to understand. Moreover, the title implies that the metasurface under consideration may achieve a multiphysics functionality. However, this does not correspond to the actual content of the paper. There is no \"multiphysics metasurface\". Instead, the authors try to determine the S-parameters of a metasurface and the temperature due to heat dissipated in the lossy substrate (but not in the conductors that are also lossy and more prone to heating).\nAs far as I understand, the paper is based on the following claims:\n1. Frequency transfer can be accomplished with the proposed combination of networks.\n2. Multiphysics analysis of metasurfaces is important.\n3. The proposed framework can be used for inverse design.\nNone of these has been shown to a degree that would merit publication. In particular:\n1. No execution time comparisons or rigorous accuracy metrics have been provided for the frequency transfer. it is unclear how to choose the bounds and the width of the extrapolated frequency band. Also, it seems that the authors' approach is data driven, although frequency transfer has recently been demonstrated for neural operators in a compelling way.\n2. The few comparisons of S_E and S_M show that there is little difference between the two. In fact, Fig. 6 shows a very low S21. Some difference between the two is observed at 6.7 GHz, but it is not explained. Is it real, or is it a numerical artifact ?\n3. Various ANNs, GANs etc have been already shown to drive inverse design of metasurfaces. Therefore, the results presented here do not demonstrate any breakthrough.", "reviewer_2": "The authors proposed a multi-physics deep learning framework that combines several techniques. This work is interesting and seems beneficial. However, the framework is complex, and each subnetwork is not novel. The manuscript does not clearly explain or compare the necessity and advantages of the different components. Moreover, the problem statement is not very clear at the beginning. To improve the paper, the reviewer suggests the authors to revise the paper by addressing the following concerns:", "reviewer_3": "In this manuscript, the authors present a data-analytic, hybrid framework that could offer zero-shot out-of-distribution (OOD) generalization with respect to operating frequency, dubbed “frequency transfer” by them. At the core of the framework is the multifidelity DeepOnet approach, which aims to incorporate both a data-driven machine learning (ML) model and a physics-based analytic model. The proposed framework is validated for both forward modeling and inverse design of electromagnetic-thermal metasurface subject to multiphysics coupling. Experimental validation shows good agreement with the computational result.\nThis manuscript addresses zero-shot out-of-distribution (OOD) generalization with respect to operating conditions, a topic that has been largely under-explored in the literature but could have far-reaching impact in wave-based systems beyond electromagnetic metasurfaces. In this regard, this manuscript could be a timely, important addition to the literature. The authors’ approach, which attempts to take advantage of both data-driven modeling and physics through the multifidelity modeling, is a good demonstration of scientific machine learning (SciML). While this point is commendable, the manuscript in the present form has major deficiencies in several different regards. It will require major revisions to meet the quality that readership of Nature Communications would anticipate.\nMajor Comments\nComment 1 (OOD generalization)\nTo raise the impact of the work, Introduction could be significantly expanded to include major approaches on OOD generalization existing in the literature. In recent years, OOD generalization has been one of the most intensely studied topics by diverse communities. Within the scope of Scientific Machine Learning (SciML), there are two mainstream approaches to extrapolation of data-driven models:\n1. Multifidelity modeling with a sparse amount of data given in the extrapolation\n2. Fine-tuning for new data\nThe two approaches in common is that they require data in the extrapolation region or re-training of models. The following reference has a good introduction of both approaches, with a good amount of empirical results and analyses included:\n- Zhu, Min, et al. \"Reliable extrapolation of deep neural operators informed by physics or sparse observations.\" Computer Methods in Applied Mechanics and Engineering 412 (2023): 116064.\nIn addition to the extrapolation approaches reported from the SciML communities, there are physics-based approaches to extrapolation, e.g., spectral generalization by wavelength normalization reported in the following reference:\n- Tanriover, Ibrahim, et al. \"Deep generative modeling and inverse design of manufacturable free-form dielectric metasurfaces.\" ACS Photonics 10.4 (2022): 875-883.\nThe authors are also recommended to go beyond the scope of frequency transfer into areas such as polarization transfer and dispersion transfer in literature review, which can help readers better understand the context of this manuscript and the relationships of these topic areas.\nComment 2 (DeepOnet)\nThe authors should consider to significantly to expand the introduction of existing methods in SciML. Most importantly, it is difficult for the reviewer to find a convincing answer to “Why DeepOnet?”. Including convincing justifications, at least qualitatively, is essential, in particular since in this manuscript no empirical validation is provided against competing techniques, e.g., physics-informed neural nets (PINN), neural operator, etc. The authors may refer to the following references which articulated in what aspects operator learning could be powerful, e.g., sample/parameter efficiency, discretization invariance, capability to capture non-local interactions:\n- Augenstein, Yannick, Taavi Repan, and Carsten Rockstuhl. \"Neural operator-based surrogate solver for free-form electromagnetic inverse design.\" ACS Photonics 10.5 (2023): 1547-1557.\n- Lee, Doksoo, et al. \"Deep Neural Operator Enabled Concurrent Multitask Design for Multifunctional Metamaterials Under Heterogeneous Fields.\" Advanced Optical Materials (2024): 2303087.\nThe authors should clearly state the exclusive advantages of including DeepOnet in the proposed framework to achieve frequency transfer.\nIn addition, considering the readership in Nature Communications, it is highly recommended to provide (1) a brief introduction to operator learning, (2) a canonical, simplified formulation on DeepOnet (as one form of supervised learning), and (3) an exposition on how it distinguishes again “conventional” machine learning surrogates, e.g., multilayer perceptron or convolutional neural nets. Such introduction centered on SciML (e.g., PINN; neural operator) seems more necessary than the reviews on VAE and GAN in Introduction.\nComment 3 (data-analytic multifidelity modeling)\nAs mentioned above, the use of multifidelity modeling as a means to incorporate machine learning and physics is commendable. Nevertheless, the descriptions and justifications on multifidelity modeling in the manuscript have a large room for improvements. The radar plot at the bottom right of Fig. 1 is better suited as a conceptual illustration of the motivation; but it currently lacks supporting arguments of the effectiveness of the approach. As a relevant point, the authors are claiming:\n“In our scheme, the pressure of multi-physics data collection is relieved by a priori knowledge, and the properties at frequency out of training range are predicted with a latent space, benefiting vast applications involving multi-frequency and multi-physics domain.”\nThe reviewer found the claim to be very weak, in particular, no quantitative investigation in any aspect of multifidelity modeling has been reported in this paper. Possible numerical experiments may include (1) performance gain against single-fidelity modeling, (2) the influence of the sampling interval ratio between Net1 and Net2, (3) the trade-off between computational cost and accuracy. The authors are encouraged to search more relevant papers built on multifidelity machine learning and include similar types of validation experiments in their manuscript. A non-exhaustive list is provided as follows:\n- Pilania, Ghanshyam, James E. Gubernatis, and Turab Lookman. \"Multi-fidelity machine learning models for accurate bandgap predictions of solids.\" Computational Materials Science 129 (2017): 156-163.\n- Howard, Amanda A., et al. \"Multifidelity deep operator networks for data-driven and physics-informed problems.\" Journal of Computational Physics 493 (2023): 112462.\n- Tran, Anh, et al. \"Multi-fidelity machine-learning with uncertainty quantification and Bayesian optimization for materials design: Application to ternary random alloys.\" The Journal of Chemical Physics 153.7 (2020).\n- Buterez, David, et al. \"Transfer learning with graph neural networks for improved molecular property prediction in the multi-fidelity setting.\" Nature Communications 15.1 (2024): 1517.\n- Li, Shibo, et al. \"Multi-fidelity Bayesian optimization via deep neural networks.\" Advances in Neural Information Processing Systems 33 (2020): 8521-8531.\n- Lyu, Yanfang, et al. \"Multi-fidelity prediction of fluid flow based on transfer learning using Fourier neural operator.\" Physics of Fluids 35.7 (2023).\nComment 4. (On the architectural design)\nIt is questionable whether the presented architectural design, illustrated in Fig. 2, consists of essential components. For example, the reviewer wonders why the same geometric input (G) should go over the progressive feature extraction twice, once in Net1 and once again in Net 2. More detailed description of the presented architecture is recommended. In doing so, naming each major net, simply indexed as Net1, Net2, … for now, in relation to their role(s) would be a good starting point.\nComment 5. (the scope of OOD generalization)\nAs opposed to in-distribution generalization that most readers would be familiar with, OOD generalization performance could involve a different scope, thus it is important to clearly specify the scope of the generalization. Throughout the manuscript, only the generalization result for the band 10-14 GHz has been reported. What would be the performance at the lower frequency regime (i.e., for a band < 2 GHz) could be included. Also important to observe is sharing whether the frequency transfer works well even beyond 14 GHz. If not, what would the deviation pattern look like? Intuitively, the error trend would presumably show a strong positive correlation w.r.t. the amount of frequency deviation (i.e., the larger the frequency deviation the larger the error would be). This could better present the strength/weakness of the presented frequency transfer.\nComment 6. (Training description)\nRegarding the architectural design and their implementation in this work, it is recommended to report the details of the architectures (e.g., the number of trainable parameters), as well as an example result of the network training (e.g., hyperparameters; deep learning platforms used for the training; optimizer) in Supporting Information. The description on the training procedure under “Numerical simulation and training details” is insufficient to reproduce what’s presented in the manuscript. While enabling reproducibility seems optional, details of the implementations could be practically informative/insightful for readers who are interested in either implementing the presented framework or developing their own upon the demonstration showcased by this work.\nComment 7. (The meta-atom Library)\nOne of the primary motivations of deep learning assisted photonics has been to embrace high-dimensional, (quasi-)freeform meta-atoms. Although the authors seem to have achieved the frequency transfer within the specified scope, the coverage of meta-atom geometry seems very narrow. Based on the reviewer’s understanding the pre-defined meta-atom representation only spans the space formed by a predefined topology, fully specified by a handful number of parameters without representation learning. Spanning and exploring only this narrow scope are questionable for the reviewer since operator learning, including DeepOnet, is known to be powerful in dealing with high-dimensional input/output fields. The authors might want to include any discussions on this point in their manuscript.\nComment 8. (Editorial proof)\nEnglish could be improved. The manuscript has a lot of typos that are not difficult to spot. To name a few:\n- Experimental experiment\n- “Take metasurface as an example” (in page 2)\n- “However, the multi-physic field coupling problems 25-27 widen the gaps between the practical engineering applications and numerical simulation.” (in page 2)\n- “Experimental experiment is also conducted to verify the accuracy of the framework with a comparison experiment to show a significant EM-thermal effects.” (in page 3)\nMore thorough proofreading is highly recommended.\nMinor Comments\nComment 9. (Descriptions on GAN)\nRegarding the comment:\n“Generative adversarial network (GAN) 39,40 , as one of the most popular frames in inverse design these years, constructs both forward and inverse models at the same time by adversarial strategy, partly alleviating the problem of non-uniqueness.”\nThe inverse model is constructed only for conditional GAN, not the vanilla GAN (Ref. 40). Vanilla GAN involves unsupervised learning (i.e., being agnostic to the property / output), thus what it does at most is “pre-processing method to compress the data space”, borrowing the authors’ expression. It has no direct contributions to both forward modeling (supervised learning) and inverse design. Thus the statement, as is, could mislead the readers who are not familiar with GAN. It is suggested to make a clear distinction between vanilla GAN, and conditional GAN as an extension.\nComment 10. (Misused terms in training)\n“is directed by three constraints”\nAll these are loss functions (the main objective of an optimization), not constraints (requirement / regularization). As another statement with a similar error:\n“Furthermore, in the area of inverse design, our approach imposes hybrid constraint to regularize the results, where the analytical model extracts useful a prior from the data with equations, and forward discriminator makes use of knowledge by build DNN, providing the paradigm of the same type of problem.”\nHerein, constraints and regularization should be rigorously distinguished.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Reply to the Reviewer 1\nDear Reviewer,\nWe sincerely thank you for providing a detailed review of our paper. We provide our replies to\nyour comments below and highlight the changes we have made in the manuscript using underlined\ntexts.\nComment 1: The paper is not well written and hence, it is difficult to understand. Moreover, the\ntitle implies that the metasurface under consideration may achieve a multiphysics functionality.\nHowever, this does not correspond to the actual content of the paper. There is no \"multiphysics\nmetasurface\". Instead, the authors try to determine the S-parameters of a metasurface and the\ntemperature due to heat dissipated in the lossy substrate (but not in the conductors that are also\nlossy and more prone to heating).\nResponse: We are sorry that we may have caused misunderstanding to you by our expressions.\nWe give our explanations and revisions from the following three aspects:\n(1) The temperature varying parameters and heat production of both substrate and conductors are\nconsidered in this work. More specifically, the relative dielectric constant of lossy substrate is\ndescried by Debye model, and the relationships between temperature and the conductivities\n(including both electrical conductivity 𝜎 and thermal conductivity 𝜅 ) are obtained from\nCu Cu\nthe experimental data in our previous publication Ref. 62. More specifically, the analytic\nfunctions of electrical conductivity (𝜎 ) and thermal conductivity 𝜅 in Cu against\nCu Cu\ntemperature is given by\n𝜎 = ∑4 𝑏 𝑇𝑖, 𝑇 ∈ (200𝐾,900𝐾) (R1)\nCu 𝑖=0 𝑖\n𝜅 = ∑1 𝑐 𝑇𝑖, 𝑇 ∈ (200𝐾,1200𝐾) (R2)\nCu 𝑖=0 𝑖\nwhere b = 2.91115×108, b = -1.5643×106, b = 3.7×103 b = -3.9347, b = 1.5644×10-3,\n0 1 2 , 3 4\nc = 420.33208, c = -0.06809. We have added these details in the revised manuscript to avoid\n0 1\nsimilar misunderstanding (Highlighted on Page 15).\nresponses of through-silicon vias applied with a periodic voltage pulse. IEEE Trans. Electron Devices 57, 1382-\n1389 (2010).\n(2) In previous version, we consider “multiphysics” in forward solver, i.e., determining the\nelectromagnetic performance (S-parameters) considering electromagnetic-thermal multi-\nphysics effects. In the revised manuscript, we further consider “multiphysics functionality” in\nthe inverse solver, i.e., optimizing metasurface structures considering both electromagnetic\nand thermal multiphysics functionality (Highlighted on Page 11).\nMore specifically, by adding temperature information (Input-T), which are the value and\nposition of maximum temperature within 2-8 GHz, to the input of inverse model, we are\ntargeted to constrain the thermal performance in metasurface besides electromagnetic\nperformance under electromagnetic-thermal coupling effects. With the same desired EM\nresponse, the achieved performances by the results from inverse model trained without Input-\nT (Option 1) and with Input-T (Option 2) are given in Fig. R1, respectively.\nFig. R1. The comparison of EM and thermal performance between the design results from inverse model\ntrained without Input-T (Option 1) and with Input-T (Option 2) under the same desired EM response. a\nand b Example I. c and d Example II.\nAs shown in Fig. R1 a and c, Option 1 and 2 have similar passband frequency\ncharacteristics, which are very close to the desired one. However, as shown in Fig. R1 b and d,\nby constraining the position of the maximum temperature near 4 GHz, the Option 2 performs\nbetter in the average temperature compared to Option 1, where the maximum temperature rise\noccurs at a higher frequency in Option 1 meaning more serious thermal effects happen.\n(3) To further verify the thermal performance by inverse method experimentally, design results in\nExample II are fabricated with 4*2 units, whose total size are 40 mm * 20 mm. As shown in\nFig. R2. a, the fabricated samples are irradiated in the waveguide with PA enhancing the\nmicrowave from VNA. After 3 mins continuous wave irradiation at 7.9-8 GHz and 20 W, the\nsurface temperature distributions of samples are given in Fig. R2. b and c, where the highest\nand average temperature of design option 1 (347.45 K, 321.22 K) are higher than those of\ndesign option 2 (331.85 K, 311.71 K), while the environment temperature is close. In addition,\nthe measured temperature at the quarter and the half of the samples (Line1, Line2) are\ncompared as shown in Fig. R2 d and e, where temperature of Option 1 (design without Input-\nT) is generally above that of Option 2 (design with Input-T).\nFig. R2. The temperature measurement of fabricated sample in waveguide. a Measurement setup. b\nInfrared image of design option 1 under 3 mins radiation. c Infrared image of design option 2 under 3\nmins radiation. d The measured temperature at the quarter of the samples. e The measured temperature at the\nhalf of the samples.\nTo conclude, in the data collection, the EM-thermal behavior of both substrate and\nconductors in the metasurface are studied in this work. Also, multi-physics functionality has\nbeen achieved by the inverse method and been verified experimentally in the revised\nmanuscript, where the temperature on metasurface at high power is reduced when EM\nperformance is satisfied. To relieve your concerns, we have added the above discussion in the\nrevised manuscript. In addition, we have also asked native speakers to help check the language\nproblem to further improve the manuscript, where we have also re-organized the literature\nreview in the Introduction Section in a more compact way to improve the organization structure.\nComment 3: As far as I understand, the paper is based on the following claims:\n1. Frequency transfer can be accomplished with the proposed combination of networks.\n2. Multiphysics analysis of metasurfaces is important.\n3. The proposed framework can be used for inverse design.\nNone of these has been shown to a degree that would merit publication.\nResponse: We have made a point-to-point reply and revision to all your concerns in the following\nresponses.\nComment 4: No execution time comparisons or rigorous accuracy metrics have been provided for\nthe frequency transfer. it is unclear how to choose the bounds and the width of the extrapolated\nfrequency band.\nResponse: We give our explanations and revisions as follows:\n(1) To quantify the overall accuracy for predicting EM response and temperature, relative error is\nused, i.e.,\n2\n𝐸 = ∑(𝑦 −𝑓(𝑥 )) ⁄∑𝑦 2 (R3)\n𝑅 𝑖 𝑖 𝑖\nwhere 𝑦 represents the true value and 𝑓(𝑥 ) is the predicted value. In predicting the average\n𝑖 𝑖\ntemperature by latent dynamics networks, 𝐸 within 2 to 10 GHz is 0.007, 𝐸 from 10.2 to 14\n𝑅 𝑅\nGHz is 0.043, which is the accuracy for the frequency transfer.\nTo quantify the prediction error at each frequency point, mean square error (MSE) is used.\nAs shown in Fig. R3 b, part of the structure samples is used for training, and other samples\nunseen in training data are used for tests, which means out-of-range tests are carried out in\nterms of both geometry and frequency ranges (frequency transfer). Compared with the results\nin Fig. R3 a, which are from traditional DNN method without the proposed Euler latent space,\nthe proposed method obtains a smaller test error, especially for out-of-range tests.\nFig. R3. The error in temperature prediction, where results above 10 GHz are frequency transfer. a\nThe results from a traditional DNN method without latent space, where 𝑓, 𝑆 , and 𝐺 are the inputs and\n𝐸\n𝑇(𝑓) is the output. b The results from the proposed latent network, where part of structure samples is used\nfor training.\n(2) For the bounds of extrapolated frequency band, on the one hand, considering the grating lobe\neffect, the application frequency of the periodic array is limited by its period, i.e., 𝑓 = 𝑐⁄2𝑑=\n𝑇\n𝑇\n𝑐 represents the speed of light, 𝑑 represents the period of array. On the other hand, as shown\nin Fig. R4a, unexpected higher order resonances occur at the frequency above 14 GHz.\nDeviations between the predicted temperature and the full wave simulation occur at these\nresonance points, as shown in Fig. R4 b. It shows the limitation of frequency transfer based on\nthe proposed method that prediction of different physical phenomena is not accurate, which\ncould be improved by combing our method with other techniques, such as transfer learning.\nHowever, due to grating lobe effect and high order resonances, there is no need to take such\nwide frequency band under consideration.\nFig. R4. Limitation of frequency transfer based on the proposed method. a Higher order resonances occur\nat the frequency above 14 GHz. b The predicted temperature compared with the results by full wave\nsimulation, which deviates at position of higher order resonances.\n(3) Further, besides the parameterized metasurface structures, we have also tested the frequency\ntransfer of proposed method on free-form meta-atoms with higher degree of freedom. More\nspecifically, we keep the framework of the prosed method unchanged but change full\nconnected network to convolutional neural network to deal with 2D patterns. As shown in Fig.\nR5, the extrapolated frequency band is set to be 11.2 to 14 GHz with 𝐸 of 0.047. It is shown\n𝑅\nthat the frequency transfers of the proposed method still validate in higher dimension cases.\nFig. R5. Average temperature by latent dynamics networks, where the meta-atoms are free-form and high-\ndimensional and the extrapolated frequency band is 11.2 to 14 GHz.\n(4) For the comparison of execution time, the EM-thermal simulation of one single geometry cost\n10 mins, considering a steady-state frequency domain solver with alternating iterations. On\ncontrary, the proposed learning-based solvers predict the results within 1 second. Further, we\nhave also included the time of accumulating dataset and training the networks in the revised\nmanuscript, where are about 7200 mins and 62 mins, respectively. It is worthwhile to mention\nthat both accumulating dataset and training the networks can be done offline and accumulating\ndataset can be further done in parallel. Further, the proposed method is also able to predict out-\nof-distributions performance with higher frequencies.\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 7) and Supplementary Information (Highlighted on Page 8-9).\nComment 5: Also, it seems that the authors' approach is data driven, although frequency transfer\nhas recently been demonstrated for neural operators in a compelling way.\nResponse: We give our explanations and revisions from the following three aspects:\n(1) Frequency transfer, which is to predict the unknown spectrum given part of the frequency data,\ncould be divided into two types according to frequency range provided by the training dataset.\na) In the first type, where the range of training dataset already cover the predicted\nfrequency: For example, the mentioned Fourier neural operator (FNO) is used to build the\nmapping relationship between the low-frequency and high-frequency wavefield in\nGeophysics fields (Ref. 36), where high-frequency data as the output is included training\ndata. Furthermore. in the field of metasurface, Ref. 37 introduces a data driven spectra-to-\nspectra strategy to decode unknown spectrum of metasurface from the latent space with the\nknown spectrum, in which both known and unknown spectrum are already used for training.\nIn Ref. 38, by combining FNO and DeepONet, the interpolation of frequency is achieved\nwith an extended FNO.\nb) The second type of frequency transfer is out-of-distribution (OOD) generalization, i.e.,\nthe predicted spectrum is not seen in the training stage, which is the problem discussed\nin this work. There are two mainstream approaches, i.e., multi-fidelity modeling and fine-\ntuning (Ref. 39). For example, transfer learning is adopted in fine-tuning FNO backbone to\nachieve OOD with a sparse amount of data given in the extrapolation (Ref. 40). In addition,\nthere are also physics-based approaches to extrapolation. In Ref. 35, based on the correlation\nof metasurface size scaling and resonant wavelength, the generalizability in operation\nwavelength in forward method is achieved. In Ref. 41, based on a prior provided by the\nequivalent circuit, spectral extrapolation of metamaterials is achieved.\nfree-form dielectric metasurfaces. ACS Photonics 10, 875-883 (2023).\nGeophys. Eng. 19, 269-282 (2022).\nelimination framework. Nat. Commun. 14, 4872 (2023).\ndata in frequency domain using neural operator. IEEE Trans. Geosci. Remote Sens. 60, 2007912 (2022).\nComputer Methods in Applied Mechanics and Engineering 412, 116064(2023).\nparametric frequency-domain wave equations. IEEE Trans. Geosci. Remote Sens. 62, 5923211 (2024).\nempowered by circuit-physics-driven deep learning. Laser & Photonics Reviews, Early Access (2024).\n(2) Despite the convincing achievement in the above literatures, the mentioned physics-based\nworks only focus on single physic response. In this work, however, the frequency transfer is\nachieved to describe a multi-physical effect, where the independent variable frequency is an\nEM quantity and the dependent variable temperature is a quantity of heat. It is noted that the\nEM and thermal behaviors in the metasurface are governed by different PDEs, and these\nPDEs are alternately solved in the steady-state solution of the EM-thermal coupling, which is hard to\ndecouple. In this case, data driven method seems to be an appropriate approach to predict EM-\nthermal behaviors.\n(3) Further, compared with the above data-driven approaches, data in the extrapolated frequency\nband are not used for training in the proposed method. Instead, the proposed method is\nsuccessfully generalized to unseen frequency bands by dynamically utilizing an Euler\nregularization in latent space.\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 2).\nComment 6: The few comparisons of S_E and S_M show that there is little difference between\nthe two. In fact, Fig. 6 shows a very low S21. Some difference between the two is observed at 6.7\nGHz, but it is not explained. Is it real, or is it a numerical artifact?\nResponse: Thank you for constructive comments. To relieve your concerns, we have fabricated\ntwo samples, one is “design sample” from the proposed method and the other is “contrast sample”.\nCompared with previous version, we have also improved our test conditions with a higher\nirradiation power and frequency to better measure the performances. We give our explanations and\nrevisions as follows:\n(1) After 10 mins irradiation under 20 W microwave of 6-7.05 GHz in waveguide, the surface\ntemperature distributions of contrast sample and design sample are given in Fig. R6. a and b,\nwhere the highest and average temperature of contrast sample (360.15 K, 342.17 K) are higher\nthan those of design sample (314.15 K, 305.12 K).\nFig. R6. Infrared image of a contrast sample and b design sample under 10 mins radiation of 6-7.05 GHz,\n(2) The increase of temperature further leads to the change of the EM properties of the sample.\nThe measured S21 recorded during irradiation are given in Fig. R7. On the one hand, as shown\nin Fig. R7 a, it is obvious that S21 of contrast sample is shifted to the left as the irradiation\ntime increase, and stable at about 9 min, where S21 at 0 min is S and that at 9 mins is S .\nE M\nOn the other hand, because the heat effect in the design sample is considered, the measured S21\nare nearly unchanged during irradiation.\nFig. R7. The change of measured S during irradiation. a Contrast sample.\n21\nb Design sample.\nTo conclude, based on the above measurement of contrast sample, it is suggested that the\nthermal effect and its changes to EM response in metasurface are significant, and it is not a\nnumerical artifact. Furthermore, if the irradiation power continues to improve, the EM-thermal\ncoupling could be much more severe and it will have a huge impact on the system stability. In\nthis way, multi-physics analysis of metasurface is important. To relieve your concerns, we\nhave added the above discussion and details in the revised manuscript (Highlighted on Page\n13).\nComment 7: Various ANNs, GANs etc have been already shown to drive inverse design of\nmetasurfaces. Therefore, the results presented here do not demonstrate any breakthrough.\nResponse: Thank you for the comment. Although, tremendous success has been witnessed in the\ninverse design of metasurface, we have made two original contributions in the inverse-design part\nof this work:\n(1) As is detailed in the Response to your Comment 2, multi-physics functionalities have been\nachieved by the inverse solver, which has been verified experimentally.\n(2) In the inverse method, since designing metasurface considering both EM and thermal\nfunctionalities is a challenging multi-objective problem, hybrid regularizations are introduced\nin inverse network to constrain and boost the training process. More specifically, besides the\nloss defined on geometry parameters, forward solvers including EM and thermal predictor and an\nanalytical model are further introduced to alleviate non-unique problem and constrain\nresonant points, respectively.\nMore specifically, we also analyze the impact of the weights of these regularizations on the\ninverse design. As shown in Fig. R8 a and b, the relative errors of these regularizations in\nvalidation set are given under different training case, in which the numbers from left to right\nin parentheses represent the weights of data-driven, analytical-driven, EM forward model, and\ntemperature forward model regularizations, respectively. E1 represents the relative error on\ngeometry determined by data-driven regularization. E2 represents the relative error on the\nresonant frequencies determined by analytical-driven regularization, in which E21 and E22\nrepresent the relative error on the first and second resonant frequency, respectively. E3\nrepresents the relative error on EM response predicted by EM forward model. E4 represents\nthe relative error on average temperature predicted by thermal forward model. Based on the\nperformances on these errors, Case5 is thought to be the most appropriate choice for our inverse\ndesign problem of multi-physic metasurface.\nFig. R8. The impact of the weights of regularizations on the inverse design. a The relative errors of\nregularizations in validation set under different training case. b Partial enlarged drawing of a. c-d The effect\nof analytical-driven regularization in this work.\nIn addition, the effect of analytical-driven regularization is visualized in Fig. R8 c and d.\nAs shown in Fig. R8 c, although EM response in the passband obtained by Case2 (governed\nby analytical-driven regularization) is less accurate than that by Case3 (governed by EM\nforward model), a higher accuracy in pole prediction is witnessed in Case2. From Fig. R8 d,\nby adding analytical model to the traditional GAN, the inverse design achieves more accurate\nlocal performance deign while ensuring the overall consistency.\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 10).\nWe hope that the explanations and revisions provided above are acceptable to you. We sincerely\nthank you once again for your valuable time and appreciate your feedback very much.\nWith thanks and best regards,\nAuthors\nReply to the Reviewer 2\nDear Reviewer,\nWe sincerely thank you for providing a detailed review of our paper. We greatly appreciate the\ntime and insight you have invested in our paper. We provide our replies to your comments below\nand highlight the changes we have made in the manuscript using underlined texts.\nGeneral Comment: The authors proposed a multi-physics deep learning framework that combines\nseveral techniques. This work is interesting and seems beneficial. However, the framework is\ncomplex, and each subnetwork is not novel. The manuscript does not clearly explain or compare\nthe necessity and advantages of the different components. Moreover, the problem statement is not\nvery clear at the beginning. To improve the paper, the reviewer suggests the authors to revise the\npaper by addressing the following concerns:\nResponse: Thank you for careful reading and constructive suggestions. We have revised the\nmanuscript to include a more detailed introduction. Also, we have conducted a complete literature\nreview and revised the introduction to clearly define the problem we are addressing. To be specific,\nwe have made a point-by-point reply and revision to all your comments below.\nComment 1: Fig.1 plots both the forward and inverse methods. Does the manuscript aim to learn\nboth? Specifically, for the forward network, G and f are inputs and T(f) are output, and for the\ninverse network, R is input and G is output. Please clearly explain the problem this work addresses.\nIn the forward method, the Multiphysics predictor is not used, and in the inverse method, the latent\ndynamics networks are not used. The reviewer suggests to rethink how to demonstrate the two\nmethods in Fig.1.\nResponse: Thank you for your constructive comments. We give our explanations and revisions as\nfollows:\nFig. R9. Procedure of MDLF, which is divided into forward method and inverse method. R represents\nrequirement, Input-T represents temperature information of metasurface, 𝐺 is geometry, 𝐺’ is the predicted\n𝐺, 𝑇 is temperature, 𝑓 is frequency, 𝑆 is S-parameter by single EM simulation, 𝑆 represents S-parameter\n𝐸 𝑀\nby Multiphysics simulation, 𝑆 ′ is the predicted 𝑆 , 𝑎(𝑓) is hidden space, 𝑓 is resonant frequency, 𝑓 ′ is\n𝑀 𝑀 𝑅 𝑅\nthe predicted 𝑓 , and 𝑆 = 𝑆 – 𝑆 .\n𝑅 𝑅 𝑀 𝐸\n(1) As presented in Fig. R9, the forward networks consist of two parts, i.e., multi-fidelity DeepOnet\nto predict electromagnetic S-parameter performance (𝑆 ) and a latent-dynamic network to\n𝑀\npredict temperature performance (𝑇(𝑓)), where both of them have considered electromagnetic\n(EM)-thermal multi-physics coupling effects.\nMore specifically, since the computational cost is high when accumulating training data\nfor predicting 𝑆 , the multi-fidelity DeepOnet is used to reduce the need of high-fidelity data\n𝑀\n(𝑆 ), where data of S-parameter performance without EM-thermal coupling effects (𝑆 ) are\n𝑀 𝐸\nfirstly accumulated to train a low-fidelity network, i.e., Net1 in Fig. R9. The output of the low-\nfidelity network is further input into the high-fidelity network (Net2) and formed as a residual\npart to reduce the difficulty in training Net2. To sum up, given geometry parameter 𝐺) and\nfrequency (𝑓) of a metasurface, the trained multi-fidelity DeepOnet can predict 𝑆 considering\n𝑀\nmulti-physics effects with reduced requirements of high-fidelity data. To validate the effects,\ncomparisons are added in the revised manuscript (Highlighted on the last two paragraphs in\n“EM response by multi-fidelity DeepONet” in Page 8).\nFurther, since obtaining thermal effects at high frequency is either computational or\nexperimental challenging, we further developed a latent-dynamic network (Net3 & Net4) to\npredict temperature performance (𝑇(𝑓)) considering EM-thermal multi-physics coupling\neffects. By using a Euler-method based dynamic latent, the proposed network can predict\ntemperature performance at higher frequency not seen in training data, i.e., out-of-distribution\n(OOD) prediction. To validate the effects of frequency transfer, results are presented in the\nrevised manuscript (Highlighted on the last paragraph in “Average temperature by latent\ndynamics networks” in Page 9).\n(2) In the inverse method, since designing metasurface considering both EM and thermal\nfunctionalities is a challenging multi-objective problem, hybrid regularizations are introduced\nin inverse network (Net5) to constrain and boost the training process. More specifically,\nbesides the loss defined on geometry parameters, forward solvers including EM and thermal\npredictor and an analytical model are further introduced to alleviate non-unique problem and\nconstrain resonant points, respectively. To validate multi-functional design and the effects of\ndifferent constrains, results are presented in the revised manuscript (Highlighted on the last\nparagraphs in Page 11).\nTo sum up, multi-physics predictor is used in the forward solvers. In previous version, thermal\nperformance is not considered in the inverse part, the strategy of latent dynamics is not used in the\ninverse solvers. In the revised manuscript, we further use both EM and thermal predictor as\nforward solver to alleviate the non-unique problem in the inverse part, which means latent\ndynamics networks is used. According to your constructive comments, we have also revised the\nexpressions and re-organized Fig.1 to avoid similar confusions in the revised manuscript.\nComment 2: In Fig.1, all networks are shown as trapezoidal boxes, which are usually used to\nrepresent encoder network. Why did the authors choose this representation? Are all networks\ndesigned with gradually decreasing dimensions?\nResponse: Thank you for careful reading and comments. In the revised manuscript, we have\nremoved these trapezoidal boxes with rectangle boxes and re-organized Fig. 1, as shown in Fig.\nR9.\nComment 3: The authors should explain why Eq.2 is referred to as the latent space and how the\na(f) contributes to better results compared to directly using G, f and S_E (f) as network inputs.\nIt is better to add the comparison result in the numerical examples.\nResponse: Thank you for the comments. We give our explanations and revisions from the\nfollowing three aspects:\n(1) Why Eq.2 is referred to as the latent space:\nAs presented in latent-dynamic network of Fig. R9, 𝑎(𝑓) is computed from the parameters\nin latent space of the whole network, where no label of 𝑎(𝑓) is provided in the whole training\nprocess and the network automatically discover a compact representation of the learning task,\ni.e., latent space. More specifically, similar with the frame of encoder and decoder, Net in Fig.\n𝐸 4\nthis way, the weights of Net are updated through the backpropagation of loss between outputs\n3\nof Net and its true value. Besides, given the initial state, i.e., 𝑎(𝑓 ) = 0, the update of latent\n4 0\nspace in frequency sequence follows Euler method, which is Eq.2 (Eq.4 now) in manuscript,\n𝜕𝑎\n𝑎(𝑓 ) = 𝑎(𝑓 )+ ×𝛽 ×∆𝑓,∆𝑓 = 𝑓 −𝑓 (R4)\n𝑛+1 𝑛 𝑛+1 𝑛\n𝜕𝑓\nwhere ∆𝑓 represents the step in the frequency sequence, 𝛽 is the hyperparameter for\n𝜕𝑎\nregularization, and is changing rate of latent space, which is also the output to Net .\n3\n𝜕𝑓\nTo conclude, given no label for the latent space, Net specifies the update rules of 𝑎(𝑓) in\n3\nthe frequency sequence, Net ensures the effectiveness of 𝑎(𝑓) in predicting the temperature\n4\nat each frequency point, which together achieve frequency transfer to the unknow frequency.\n(2) How the a(f) contributes to better results: The latent dynamics network has demonstrated\nits effectiveness in Ref. 60 for predicting the evolution, even in time-extrapolation scenarios,\nof dynamic system. In this work, we introduced this strategy for frequency transfer, where\ninstead of predicting temporal sequences, frequency sequences are predicted. The effectiveness\nis due to the following reasons:\n(a) Firstly, as also pointed out by Ref. 60, the latent dynamics network enables the training\nalgorithm to select a compact representation of the state that is functional in predicting\nits dynamics (i.e., frequency sequences in our work). An auto-encoder, conversely,\nwhen trained, extracts feature on a purely statistical basis, being agnostic of the\nimportance of each feature in determining dynamic system. In other words, the latent\nstates allow propagate nonlocal information across different frequency domains.\n(b) Secondly, the proposed architecture enables the sharing of the trainable parameters\nacross different frequency points and geometry parameters, where it has been\ndemonstrated by our previous work (Ref. 55) and other literatures (Ref. 56) that weight\nsharing is extremely important for improving generalization abilities, i.e., frequency\ntransfer in this work, of networks.\napproaches by multi-domain weight-sharing for full-wave inverse scattering. IEEE Trans. Geosci.\nRemote Sens. 61, DOI:10.1109/TGRS.2023.3303572 (2023).\n2016).\nof spatio-temporal processes through latent dynamics networks. Nat. Commun. 15, 1834 (2024).\n(c) The trajectories of the latent space are also an important factor. In Fig. R10 of this reply\nletter, given the training data from 2-10 GHz, we further present the evolution of latent\nspace from the initial frequency 𝑓 = 2 GHz to the maximum observation frequency 14\ndomain, where the smooth changes also boost extrapolation predictions (Ref. 60),\ncompared with disordered latent distributions.\nFig. R10. The evolution of latent space from 2 GHz to 14 GHz, where the spectrum within 2 to 10 GHz is\ngiven in training set and the prediction from 10.2 to 14 GHz is frequency transfer, where different colors\nrepresent different frequencies.\n(3) Comparison results:\n(a) To evaluate the contribution of 𝑎(𝑓) in temperature prediction as well as frequency transfer,\nwe added a contrast test between the traditional method without the proposed latent space\nand the proposed method with latent space. In the traditional method, 𝑓, 𝑆 , and 𝐺 are the\n𝐸\ninputs to DNN and 𝑇(𝑓) is the output, which has a similar structure with Net2 but has a\ndifferent output.\nSpecifically, to quantify the overall accuracy for predicting temperature considering\nEM-thermal coupling effects, relative error (𝐸 ) is used. In proposed method, 𝐸 within 2\n𝑅 𝑅\nto 10 GHz is about 0.007, 𝐸 from 10.2 to 14 GHz is about 0.043, which is frequency\n𝑅\ntransfer. In the traditional method, with the same dataset, 𝐸 within 2 to 10 GHz is about\n𝑅\n0.0281, and that in 10.2 to 14 GHz is about 0.1643, which is larger than the error in the\nproposed method.\nFig. R11. The error in temperature prediction, where results above 10 GHz are frequency transfer. a The\nresults from a traditional DNN method without latent space, where 𝑓, 𝑆 , and 𝐺 are the inputs and\n𝐸\n𝑇(𝑓) is the output. b The results from the proposed latent dynamics network, where part of structure samples is\nused for training.\nTo quantify the prediction error at each frequency point, mean square error (MSE) is\nused. As shown in Fig. R11 b, in the proposed model, part of the structure samples is used\nfor training, which means out-of-range tests are carried out in both geometry and frequency\nrange. Compared with the results in Fig. R11 a, which are from a traditional DNN method\nwithout latent space, the proposed method obtains a smaller test error.\n(b) In addition, some examples of temperature prediction are given in Fig. R12, where the\ncomparisons between the full-wave simulation, the proposed method, and the traditional\nmethod are conducted. From the figures, even if the curves by the traditional method meet\nthe trends of the simulated results, the predicted temperature deviate from the simulation,\nwhen it comes to the extrapolated frequency points.\nFig. R12 The comparison between the average temperature by simulation, that predicted by the proposed\nmethod, and that by the traditional DNN method where 𝑓, 𝑆 , and 𝐺 are the inputs and 𝑇(𝑓) is the output.\n𝐸\nTo relieve your concerns, we have added the above comparisons and discussions in the revised\nmanuscript (Highlighted on Page 5 and 9) and Supplementary Information (Highlighted on Page\n8).\nComment 4: The discrete frequencies in Fig.2(a) and Fig.2(c) are denoted as f(n) and f(i),\nrespectively. If they are the same, please use consistent notation.\nResponse: Thank you for comments. To avoid confusion, we have unified the notation in the\nrevised manuscript.\nComment 5: The authors should add the legends for the lines in Fig.3(a) and explain the meanings\nof blue and right arrows in Fig.3\nResponse: Thank you for your comments. The arrows in Fig.3 connect the different lines and their\ncorresponding vertical axis. In the revised manuscript, we have added the legends for the lines and\nthe meanings of blue and right arrows in Fig.3.\nComment 6: Fig.3(b) compares results from simulations, Net 1 with low-fidelity data and Net 1+\nNet 2 with both low-fidelity data and high-fidelity data, demonstrating the accuracy of Net1+Net2.\nHowever, the reviewer suggests the authors add the results based on the Network only with high-\nfidelity data (data generated by EM-thermal simulation) to further illustrate the importance of\nmulti-fidelity DeepONet.\nResponse: Thank you for your constructive comments. According to your suggestion, the Network\ntrained only with high-fidelity data (single-fidelity DeepONet) is compared with the proposed\nmulti-fidelity model to illustrate the importance of the proposed multi-fidelity DeepONet. More\nspecifically, both the trained single-fidelity DeepONet and multi-fidelity DeepONet are tested with\n106 random samples that are unseen in the training process, where the multi-fidelity DeepONet\nobtains an average relative error of 0.036, while the relative error of single-fidelity DeepONet is\n0.154. Therefore, with a preliminary knowledge of space mapping provided by low-fidelity model\nNet1, the multi-fidelity strategy allows a relatively smaller high-fidelity dataset, which is also\nhelpful in other complex scenarios of field coupling. In the revised manuscript, we have added the\nquantitative comparisons according to your constructive suggestions (Highlighted on Page 8).\nComment 7: Fig.3(b) shows the accuracy of simulation and the proposed networks. However, the\nmanuscript does not discuss the advantages of network over simulation. If simulation is time-\nconsuming, include a comparison of computational time.\nResponse: Thank you for the comment. In the revised manuscript, we have added quantitative\ncomparisons (Highlighted on Page 16). More specifically, the EM-thermal simulation of one single\ngeometry cost 10 mins, considering a steady-state frequency domain solver with alternating\niterations. On contrary, the proposed learning-based solvers predict the results within 1 second.\nFurther, we have also included the time of accumulating dataset and training the networks in the\nrevised manuscript, where are about 7200 mins and 62 mins, respectively. It is worthwhile to\nmention that both accumulating dataset and training the networks can be done offline and\naccumulating dataset can be further done in parallel. Further, the proposed method is also able to\npredict out-of-distributions performance with higher frequencies.\nComment 8: In Fig.4(c), the results for Cases 1, 3 and 5 are worse than those for Cases 2 and 4,\nwhich are similar. Cases 1, 3 and 5 share w2=0. Explain why the analytical model does not improve\naccuracy when the forward model is present. It seems the analytical model is not necessary.\nResponse: For your concern about the necessity of analytical model in our proposed method, we\ngive our revisions and explanations from the following three aspects:\n(1) Due to a small range of weight adjustment strategy, the performance in different cases was\nclose in the previously submitted manuscript. In the revised manuscript, we have adopted a\nlarger range of weights, in which the average temperature predicted by the thermal forward\nmodel is also under investigation.\nAs shown in Fig. R13 a and b, the relative errors in validation set are given under different\ntraining case, in which the numbers in Fig. R13 a from left to right in parentheses represent the\nweights of data-driven regularization w , analytical model regularization w , EM forward\n1 2\nmodel regularization w , and temperature forward model regularizations w , respectively.\n3 4\nBesides, E1 represents the relative error on geometry determined by data-driven regularization.\nE21 and E22 represent the relative error on the first and second resonant frequency determined\nby analytical-driven regularization, respectively. E3 represents the relative error on EM\nresponse predicted by EM forward model. E4 represents the relative error on average\ntemperature predicted by thermal forward model.\nFig. R13. The impact of the weights of regularizations on the inverse design. a The relative errors of\nregularizations in validation set under different training case. b Partial enlarged drawing of a.\n(2) As shown in Fig. R13 b, compared with Case3 driven by data and EM forward model, Case2\nwhich is driven by data and EM forward model, has a similar E21 and a lower E22, which\nmeans a better prediction in resonant frequency. From Case1 to Case3 and Case5, EM forward\nmodel and analytical model are added into traditional method in turns, in which a gradual\naccuracy improvement can be observed. In addition, Table R1 gives the detailed value in Fig.\nR13, where the lowest value in each evaluation index is bolded. From the table, Case5, which\nis adopted as the proposed data-analytical driven model, performs the best in many indexes\nand has very close value with the best performance in other index.\nTable R1. The relative error under cases with different weights\nCase w w w w E1 E21 E22 E3 E4\n1 2 3 4\n1 1 0 0 0 0.315 0.0624 0.0681 0.141 0.0341\n2 0.2 0.8 0 0 0.335 0.0551 0.0600 0.151 0.0327\n3 0.3 0 0.7 0 0.294 0.0547 0.0672 0.0907 0.0305\n4 0.2 0 0 0.8 0.291 0.0610 0.0680 0.140 0.0290\n5 0.15 0.5 0.35 0 0.294 0.0539 0.0597 0.0811 0.0254\n6 0.2 0.4 0 0.4 0.302 0.559 0.0587 0.112 0.0288\n7 0.1 0.4 0.2 0.3 0.310 0.0547 0.0600 0.0893 0.0257\nFig. R14. The effect of analytical-driven regularization in this work.\n(3) Furthermore, the effect of analytical-driven regularization is visualized in Fig. R14 a and b. As\nshown in Fig. R14 a, although EM response in the passband obtained by Case2 (governed\nby analytical model) is less accurate than that by Case3 (governed by EM forward model), a higher\naccuracy in pole prediction is witnessed in Case2. From Fig. R14 b, by adding analytical model\nto the model driven by data and EM forward method, the inverse design achieves more accurate\nlocal performance deign while ensuring the overall consistency.\nTo relieve your concerns, we have added the above information in the revised manuscript\n(Highlighted on Page 10).\nComment 9: Fig.4(c) shows Case 2 has a lower error of S_M than Case 4, but a slightly higher\nerror for f_R. Why was Case 4 chosen as the final configuration? Is the accuracy f_R more critical\nthan G? The authors need to clarify which quantities are most important for the design\nResponse: Thank you for careful reading. In choosing the weights in the development of inverse\ndesign network, a comprehensive consideration of all the performance factors is taken. As shown\nin Table R1, Case5, which is adopted as the proposed data-analytical driven model, performs the\nbest in most indexes and has very close value with the best performance in other index. To relieve\nyour concerns, we have added the above information in the revised manuscript and Supplementary\nInformation (Highlighted on Page 10).\nComment 10: In Fig.4b, why does the inverse-designed EM response have a peak at the end of\nthe frequency?\nResponse: We give our explanations and revisions as follows:\n(1) In the inverse design scenario, usually we do not have detailed performance distributions, i.e.,\ndesired EM performance, before we design. Instead, we usually have some incomplete\nrequirement, for example, we want to find a geometry with two passbands at about 4 GHz and\n8-10 GHz, thus a polyline of performance is draw according to these obscure requirements, as\nshown in Fig. R15a (Fig.4b in the manuscript), indicating the desired passband location and\nbandwidth. However, no one knows whether metasurface structure satisfying the plotted\npolyline exists or not. In this case, the aim of inverse method here is to find a geometry to\nsatisfy the fictitious performance as much as possible, deviations outside the required passband\ncould be acceptable. Besides, re-trained with weights of Case5 in Fig. R15, the results by Net\n5\ncould be updated as shown in Fig. R15b, in which the increase at the end of frequency has been\ndecreased but not eliminated, considering the EM properties of the geometry itself. In the\nrevised manuscript, we have replaced Fig.4b (Fig. 4e now) with Fig. R15.\nFig. R15. The comparison between the desired EM response and that inversely designed by Net , where\n5\nthe desired one is a polyline connected with a few points available. a Results by the previous Net .\n5\nb Results\nby the re-trained Net .\n5\n(2) Furthermore, in the situation where real EM response curves from simulation are available,\nconsistency between the desired EM response and the EM response from the designed\ngeometry is observed, as shown in Fig. R16. Particularly, in Fig. R16b, the characteristic of a\npeak at the end of the frequency is captured by the inverse method, and the inversely designed\ngeometry can reproduce this feature.\nFig. R16. The comparison between the desired EM response and that inversely designed by Net , where\n5\nthe desired one is obtained by simulation.\nThank you for your constructive comments, we have added the above discussions in the revised\nmanuscript (Highlighted on Page 10-11).\nComment 11: In Fig.5(c), the measured S_E is compared with simulation. Include the measured\nS_M in the same figure to highlight the importance of Net2.\nResponse: Thank you for your constructive comments. We give our explanations and revisions as\nfollows:\n(1) The numerical simulation of 𝑆 is based on frequency domain EM-thermal algorithm, which\n𝑀\nmeans the measurement of 𝑆 should include long-time, high-power irradiations to the sample\n𝑀\nat each single frequency point. Besides, the measurement should be under real-time\nmonitoring to ensure the sample reaches a stable state. Consequently, experimental conditions\ncompletely consistent with the simulation can hardly be achieved, especially for a high-power\nmicrowave source that can cover a wide band from 2 to 14 GHz enabling space radiation to\nreach the simulation level.\n(2) Instead, to verified Net2, we have compared the predicted results with full-wave multi-physical\nsimulation in Fig. 3. Further, to highlight the importance of Net2 experimentally, we have\nconducted the following two experiments:\na) The measurement of the EM response under high power irradiation is conducted in the\nwaveguide with power amplifier (PA) enhancing the microwave from VNA as shown in\nFig. R17a. The inversely designed result is fabricated with 4*2 units with the size of 40\nmm * 20 mm, where the desired target is to fulfil the EM response as well as reducing the\nheat effects. After 9 mins continuous wave irradiation at 6-7.05 GHz and 20 W, the\nmeasured S21 are basically unchanged during irradiation, since the heat effect in the\ndesign sample is designed to be NOT severe in the observed frequency, whose highest\nand average temperature after 10 mins irradiation are 314.15 K, 305.12 K. However, due\nto the difference in the EM behaviors of small samples in waveguides and periodic large\narrays in free space, the measured S21 in waveguide is not compared with that in Fig.\n5(c), which is the measurement in free space.\nFig. R17. The measurement of the EM response under high power irradiation in the waveguide. a\nMeasurement setup. b The measurement results of the designed sample.\na) To further highlight the importance of Net2, which is the importance of predicting the\nchanges of EM response under multi-physical coupling, a contrast sample is further\nfabricated with the same size of designed sample and measured in the condition of Fig.\nR17a. As shown in Fig. R18a, after 10 mins irradiation under 20 W microwave of 6-7.05\nGHz in waveguide, the surface temperature distribution of contrast sample is given, where\nthe highest and average temperature is 360.15 K and 342.17 K. The increase of\ntemperature further leads to the change of the EM properties of the sample. As shown in\nFig. R18 b, it is obvious that S21 of contrast sample is shifted to the left as the irradiation\ntime increase, and basically stable at about 9 min, where S21 at 0 min is 𝑆 and that at 9\n𝐸\nmins is 𝑆 . Therefore, it is seen that EM response is significantly affected by multi-\n𝑀\nphysical coupling, and Net2 that considers multi-physical coupling is important.\nFig. R18. The EM and thermal performance of contrast sample under high power irradiation. a Infrared\nimage under 10 mins radiation of 6-7.05 GHz, 20 W microwave in waveguide. b The change of measured\nS21 during irradiation.\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 13).\nComment 12: There is a typo on Page 8: it should be “the forward discriminator” instead of “he\nforward discriminator”.\nResponse: Thank you very much for the careful check. We have corrected the mentioned typos in\nthe revised manuscript.\nComment 13: What is S_E^' in Page 3? Is it the same as S_E in Fig.1?\nResponse: 𝑆 ′ is the predicted 𝑆 . In the revised manuscript, we have added this information in\n𝐸 𝐸\nthe illustration of Fig. 1.\nComment 14: The vanilla DeepONet by Lu et al. merges branch and trunk outputs through a dot\nproduct operation. The authors state that their output S_E^' is generated by the convolution\nbetween the two last layers. Explain the modifications made to DeepONet and the reasons behind\nthese changes.\nResponse: We are sorry to cause misunderstanding to you by our expressions. We also merge\nbranch and trunk outputs through a dot product operation. In the revised manuscript, we have\ncorrected this typo.\nComment 15: The authors should introduce key terminologies in the introduction, such as S-\nparameter, single-EM simulation, multi-physical simulation.\nResponse: Thank you for the suggestion. S-parameter is the scattering parameters, which reflects\nthe characteristics of the reflected/transmitted signal in the frequency domain. Single-EM\nsimulation refers to numerical calculation of object/system’s electromagnetic (EM) behaviors,\nwhich is governed by Maxwell’s equations. Multi-physical simulation in this work, distinguished\nfrom single-EM simulation, refers to the numerical calculation that consider both the EM and\nthermal multi-physical behaviors of an object or system, which is governed by Maxwell’s\nequations and heat conduction equation. Here, multi-physical behaviors in metasurface include\nEM propagation, heat loss, and heat transfer, which is expressed as S-parameter and temperature.\nTo relieve your concerns, we have added the above information in the introduction in the revised\nmanuscript (Highlighted on Page 2).\nWe hope that the explanations and revisions provided above are acceptable to you. We sincerely\nthank you once again for your valuable time and appreciate your feedback very much.\nWith thanks and best regards,\nAuthors\nReply to the Reviewer 3\nDear Reviewer,\nWe sincerely thank you for providing a detailed review of our paper. We greatly appreciate the\ntime and insight you have invested in our paper. We provide our replies to your comments below\nand highlight the changes we have made in the manuscript using underlined texts.\nGeneral Comment: In this manuscript, the authors present a data-analytic, hybrid framework that\ncould offer zero-shot out-of-distribution (OOD) generalization with respect to operating frequency,\ndubbed “frequency transfer” by them. At the core of the framework is the multifidelity DeepOnet\napproach, which aims to incorporate both a data-driven machine learning (ML) model and a\nphysics-based analytic model. The proposed framework is validated for both forward modeling\nand inverse design of electromagnetic-thermal metasurface subject to multiphysics coupling.\nExperimental validation shows good agreement with the computational result.\nThis manuscript addresses zero-shot out-of-distribution (OOD) generalization with respect to\noperating conditions, a topic that has been largely under-explored in the literature but could have\nfar-reaching impact in wave-based systems beyond electromagnetic metasurfaces. In this regard,\nthis manuscript could be a timely, important addition to the literature. The authors’ approach,\nwhich attempts to take advantage of both data-driven modeling and physics through the\nmultifidelity modeling, is a good demonstration of scientific machine learning (SciML). While\nthis point is commendable, the manuscript in the present form has major deficiencies in several\ndifferent regards. It will require major revisions to meet the quality that readership of Nature\nCommunications would anticipate.\nResponse: Thank you very much for your recognition of our work. We have made a point-to-point\nreply and revision to all your concerns in the following responses.\nComment 1: To raise the impact of the work, Introduction could be significantly expanded to\ninclude major approaches on OOD generalization existing in the literature. In recent years, OOD\ngeneralization has been one of the most intensely studied topics by diverse communities. Within\nthe scope of Scientific Machine Learning (SciML), there are two mainstream approaches to\nextrapolation of data-driven models:\n1. Multifidelity modeling with a sparse amount of data given in the extrapolation\n2. Fine-tuning for new data\nThe two approaches in common is that they require data in the extrapolation region or re-training\nof models. The following reference has a good introduction of both approaches, with a good\namount of empirical results and analyses included:\n- Zhu, Min, et al. \"Reliable extrapolation of deep neural operators informed by physics or sparse\nobservations.\" Computer Methods in Applied Mechanics and Engineering 412 (2023): 116064.\nIn addition to the extrapolation approaches reported from the SciML communities, there are\nphysics-based approaches to extrapolation, e.g., spectral generalization by wavelength\nnormalization reported in the following reference:\n- Tanriover, Ibrahim, et al. \"Deep generative modeling and inverse design of manufacturable free-\nform dielectric metasurfaces.\" ACS Photonics 10.4 (2022): 875-883.\nThe authors are also recommended to go beyond the scope of frequency transfer into areas such as\npolarization transfer and dispersion transfer in literature review, which can help readers better\nunderstand the context of this manuscript and the relationships of these topic areas.\nResponse: Thank you for your constructive comments. We have re-organized the Introduction to\nhelp readers better understand the context of this manuscript and cited the above-mentioned\nliterature in the revised manuscript. To be specific:\nNowadays, generalizations with respect to operating state of metasurface, like polarization\ntransfer and dispersion transfer (Ref. 35), have attracted many attentions to relieve the burden of\naccumulating dataset and improve applicability of model. Frequency transfer, which is to predict\nthe unknown spectrum given part of the frequency data, could be divided into two types according\nto frequency range provided by the training dataset.\na) In the first type, the range of training dataset already covers the predicted frequency.\nFor example, the mentioned Fourier neural operator (FNO) is used to build the mapping\nrelationship between the low-frequency and high-frequency wavefield in Geophysics fields\n(Ref. 36), where high-frequency data as the output is included training data. Furthermore. in\nthe field of metasurface, Ref. 37 introduces a data driven spectra-to-spectra strategy to\ndecode unknown spectrum of metasurface from the latent space with the known spectrum,\nin which both known and unknown spectrum are already used for training. In Ref. 38, by\ncombining FNO and DeepONet, the interpolation of frequency is achieved with an extended\nFNO.\nb) The second type of frequency transfer is out-of-distribution (OOD) generalization, i.e.,\nthe predicted spectrum is not seen in the training stage, which is the problem discussed\nin this work. There are two mainstream approaches, i.e., multi-fidelity modeling and fine-\ntuning (Ref. 39). For example, transfer learning is adopted in fine-tuning FNO backbone to\nachieve OOD with a sparse amount of data given in the extrapolation (Ref. 40). In addition,\nthere are also physics-based approaches to extrapolation. In Ref. 35, based on the correlation\nof metasurface size scaling and resonant wavelength, the generalizability in operation\nwavelength in forward method is achieved. In Ref. 41, based on a prior provided by the\nequivalent circuit, spectral extrapolation of metamaterials is achieved.\nfree-form dielectric metasurfaces. ACS Photonics 10, 875-883 (2023).\nGeophys. Eng. 19, 269-282 (2022).\nelimination framework. Nat. Commun. 14, 4872 (2023).\ndata in frequency domain using neural operator. IEEE Trans. Geosci. Remote Sens. 60, 2007912 (2022).\nComputer Methods in Applied Mechanics and Engineering 412, 116064(2023).\nparametric frequency-domain wave equations. IEEE Trans. Geosci. Remote Sens. 62, 5923211 (2024).\nempowered by circuit-physics-driven deep learning. Laser & Photonics Reviews, Early Access (2024).\nc) Despite the convincing achievement in the above literatures, the mentioned physics-based\nworks only focus on single physic response. In this work, however, the frequency transfer is\nachieved to describe a multi-physical effect, where the independent variable frequency is an\nEM quantity and the dependent variable temperature is a quantity of heat. It is noted that the\nEM and thermal behaviors in the metasurface are governed by different PDEs, and these\nPDEs are alternately solved in the steady-state solution of the EM-thermal coupling, which\nis hard to decouple. In this case, data driven method seems to be an appropriate approach to\npredict EM-thermal behaviors.\nFurther, compared with the above data-driven approaches, data in the extrapolated\nfrequency band are not used for training in the proposed method. Instead, the proposed\nmethod is successfully generalized to unseen frequency bands by dynamically utilizing Euler\nregularization in latent space.\nComment 2: The authors should consider to significantly to expand the introduction of existing\nmethods in SciML. Most importantly, it is difficult for the reviewer to find a convincing answer\nto “Why DeepOnet?”. Including convincing justifications, at least qualitatively, is essential,\nin particular since in this manuscript no empirical validation is provided against competing\ntechniques, e.g., physics-informed neural nets (PINN), neural operator, etc. The authors may refer\nto the following references which articulated in what aspects operator learning could be powerful,\ne.g., sample/parameter efficiency, discretization invariance, capability to capture non-local\ninteractions:\n- Augenstein, Yannick, Taavi Repan, and Carsten Rockstuhl. \"Neural operator-based surrogate\nsolver for free-form electromagnetic inverse design.\" ACS Photonics 10.5 (2023): 1547-1557.\n- Lee, Doksoo, et al. \"Deep Neural Operator Enabled Concurrent Multitask Design for\nMultifunctional Metamaterials Under Heterogeneous Fields.\" Advanced Optical Materials (2024):\n2303087.\nThe authors should clearly state the exclusive advantages of including DeepOnet in the proposed\nframework to achieve frequency transfer.\nIn addition, considering the readership in Nature Communications, it is highly recommended to\nprovide (1) a brief introduction to operator learning, (2) a canonical, simplified formulation on\nDeepOnet (as one form of supervised learning), and (3) an exposition on how it distinguishes again\n“conventional” machine learning surrogates, e.g., multilayer perceptron or convolutional neural\nnets. Such introduction centered on SciML (e.g., PINN; neural operator) seems more necessary\nthan the reviews on VAE and GAN in Introduction.\nResponse: Thank you for your constructive comments. We give our explanations and revisions\nabout using the mentioned network as follows:\n(1) Sample/Parameter efficiency: A DeepONet consists of two sub-networks, i.e., branch net\nand trunk net, where branch net encodes the input function at a fixed number and trunk net\nencodes the locations for the output functions (Ref. 57). Consequently, by taking frequency as\nthe input to the trunk net, DeepONet is a suitable frame to achieve a free sampling of frequency\nand a continuous representation of the spectrum in metasurface.\nIn this work, the advantages of sample/parameter efficiency of neural operators enable the\nsharing of the trainable parameters across different frequency points and geometry\nparameters, where it has been demonstrated by our previous work (Ref. 55) and other\nliteratures (Ref. 56) that weight sharing is extremely important for improving generalization\nabilities, i.e., frequency transfer in this work, of networks. To validate the effects, the effects\nof frequency transfer are presents in the revised manuscript (Fig. 3c).\n(2) Capability to capture non-local interactions: Combing DeepONet with latent dynamics in\nour work, enables the training algorithm to select a compact representation of the state that is\nfunctional in predicting its dynamics (i.e., frequency sequences in our work). In other words,\nthe latent states allow propagate nonlocal information across different frequency domains,\nwhich is important for out-of-distribution predictions, i.e., frequency transfer in this work.\n(3) Small Dataset Training: It is also reported that DeepONets could learn operators accurately\nand efficiently from a relatively small dataset (Ref. 57), which is very helpful in predicting the\nelectromagnetic (EM)-thermal effects of metasurface, since obtaining thermal effects is either\ncomputational or experimental challenging, especially when it comes to high frequencies.\nMore specifically, in this work, since the computational cost is high when accumulating\ntraining data with Multiphysics coupling effects, the multi-fidelity DeepOnet is used to reduce\nthe need of high-fidelity data, where data of S-parameter performance without EM-thermal\ncoupling effects are firstly accumulated to train a low-fidelity network. The output of the low-\nfidelity network is further input into the high-fidelity network and formed as a residual part to\nreduce the difficulty in training.\nTo validate the effects of multi-fidelity DeepOnet, the Network trained only with high-\nfidelity data (single-fidelity DeepONet) is compared with the proposed multi-fidelity model to\nillustrate the importance of the proposed multi-fidelity DeepONet, where details have been\npresented in the response to your next comments.\n(4) Complex Physical Phenomena Predictions: Nowadays, deep neural operators have shown\ntremendous success in modelling complex physical phenomena and engineering problem (Ref.\n11, 53, 54), where the data efficiency of surrogate solver are improved in learning the\nrelationship between the high-dimensional structure and response. In the revised manuscript,\nwe have further added the predictions of free-form meta-atoms with higher degree of freedom\nby the proposed method (Highlighted on the last paragraph in “Average temperature by latent\ndynamics networks” in Page 9).\nTo relieve your concerns, we have added the above information in the introduction in the\nrevised manuscript. Also, an introduction of SciML, especially deep neural operator, have been\nadded, as well as their simplified formulations and applications.\nmultifunctional metamaterials under heterogeneous fields. Adv. Opt. Mater. (2024).\ncharacterization and inverse design of stochastic architected metamaterials using neural operators.\narXiv preprint arXiv: 2311.13812 (2023).\ninverse design. ACS Photonics 10, 1547-1557. (2023).\nmulti-domain weight-sharing for full-wave inverse scattering. IEEE Trans. Geosci. Remote Sens. 61,\nDOI:10.1109/TGRS.2023.3303572 (2023).\n2016).\nequations based on the universal approximation theorem of operators. arXiv preprint arXiv:1910.03193 (2019).\ntemporal processes through latent dynamics networks. Nat. Commun. 15, 1834 (2024).\nComment 3: As mentioned above, the use of multifidelity modeling as a means to incorporate\nmachine learning and physics is commendable. Nevertheless, the descriptions and justifications\non multifidelity modeling in the manuscript have a large room for improvements. The radar plot\nat the bottom right of Fig. 1 is better suited as a conceptual illustration of the motivation; but it\ncurrently lacks supporting arguments of the effectiveness of the approach. As a relevant point, the\nauthors are claiming:\n“In our scheme, the pressure of multi-physics data collection is relieved by a priori knowledge,\nand the properties at frequency out of training range are predicted with a latent space, benefiting\nvast applications involving multi-frequency and multi-physics domain.”\nThe reviewer found the claim to be very weak, in particular, no quantitative investigation in any\naspect of multifidelity modeling has been reported in this paper. Possible numerical experiments\nmay include (1) performance gain against single-fidelity modeling, (2) the influence of the\nsampling interval ratio between Net1 and Net2, (3) the trade-off between computational cost and\naccuracy. The authors are encouraged to search more relevant papers built on multifidelity\nmachine learning and include similar types of validation experiments in their manuscript. A non-\nexhaustive list is provided as follows:\n- Pilania, Ghanshyam, James E. Gubernatis, and Turab Lookman. \"Multi-fidelity machine learning\nmodels for accurate bandgap predictions of solids.\" Computational Materials Science 129 (2017):\n156-163.\n- Howard, Amanda A., et al. \"Multifidelity deep operator networks for data-driven and physics-\ninformed problems.\" Journal of Computational Physics 493 (2023): 112462.\n- Tran, Anh, et al. \"Multi-fidelity machine-learning with uncertainty quantification and Bayesian\noptimization for materials design: Application to ternary random alloys.\" The Journal of Chemical\nPhysics 153.7 (2020).\n- Buterez, David, et al. \"Transfer learning with graph neural networks for improved molecular\nproperty prediction in the multi-fidelity setting.\" Nature Communications 15.1 (2024): 1517.\n- Li, Shibo, et al. \"Multi-fidelity Bayesian optimization via deep neural networks.\" Advances in\nNeural Information Processing Systems 33 (2020): 8521-8531.\n- Lyu, Yanfang, et al. \"Multi-fidelity prediction of fluid flow based on transfer learning using\nFourier neural operator.\" Physics of Fluids 35.7 (2023).\nResponse: Thank you for your constructive comments. We give our explanations and revisions as\nfollows:\n(1) According to your suggestion, the Network trained only with high-fidelity data (single-fidelity\nDeepONet) is compared with the proposed multi-fidelity model to illustrate the importance of\nthe proposed multi-fidelity DeepONet. More specifically, both the trained single-fidelity\nDeepONet and multi-fidelity DeepONet are tested with 106 random samples that are unseen in\nthe training process, where the multi-fidelity DeepONet obtains an average relative error of\n0.036, while the relative error of single-fidelity DeepONet is 0.154. Therefore, with a\npreliminary knowledge of space mapping provided by low-fidelity model Net1, the multi-\nfidelity strategy allows a relatively smaller high-fidelity dataset, which is also helpful in other\ncomplex scenarios of field coupling. In the revised manuscript, we have added the quantitative\ncomparisons according to your constructive suggestions.\n(2) By adjusting the sampling strategy of high-fidelity data, the influence of sampling interval ratio\nbetween Net1 and Net2 has been studied in the revised manuscript. As indicated in the left of\nTable R2, as the sampling interval ratio increases from 2:1 to 12:1, the relative error in\nvalidation set rise like a step. For example, the error in Case-R2 and Case-R3 is very close,\nhowever, from Case-R1 to Case-R2 and to Case-R4, the rise in error is more obvious, since\ntheir interval ratio increases exponentially. A similar situation can be found in Case-R3 and\nCase-R6. In addition, if the training data is randomly selected in the original high-fidelity data,\nwhich has 61 points in one single geometry, the influences of the sampling number in high-\nfidelity data are shown on the right of Table R2.\nTable R2. The influence of sampling strategy to the relative error in validation set\nInterval ratio between Relative error in Number of high- Relative error in\nNet1 and Net2 validation set fidelity data validation set\n2:1 (Case-R1) 0.036 61 0.036\n4:1 (Case-R2) 0.042 31 0.040\n6:1 (Case-R3) 0.043 21 0.040\n8:1 (Case-R4) 0.048 11 0.045\n10:1 (Case-R5) 0.050 6 0.050\n12:1 (Case-R6) 0.071 - -\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 8) and cited the above literatures to enrich literature review.\nComment 4: It is questionable whether the presented architectural design, illustrated in Fig. 2,\nconsists of essential components. For example, the reviewer wonders why the same geometric\ninput (G) should go over the progressive feature extraction twice, once in Net1 and once again in\nNet 2. More detailed description of the presented architecture is recommended. In doing so,\nnaming each major net, simply indexed as Net1, Net2, … for now, in relation to their role(s) would\nbe a good starting point.\nResponse: Thank you for the comment. We give our explanations and revisions as follows:\n(1) On the one hand, in vanilla multi-fidelity DeepONet proposed by Lu. et al. (Ref. 61), it is\nrecommended to append evaluation of low-fidelity in the trunk-net inputs and leave the branch-\nnet unchanged, i.e.,\n𝒢 (𝑣)(𝜉) = 𝒬(𝑣)(𝜉,𝒢 (𝑣)(𝜉)) (R5)\n𝐻 𝐿\nwhere 𝒢 represents the high-fidelity model, 𝒢 represents the low-fidelity model. 𝑣\n𝐻 𝐿\nrepresents the input function, which is the geometric input (𝐺) in this work. 𝜉 represents the\nlocation of the output function, which is frequency input (𝑓) in this work. Therefore, 𝑣 as the\nonly input to the branch-net is essential in both high-fidelity and low-fidelity model. With\nreference to the above framework, the geometric input (𝐺) is fed twice into the network, once\nin Net1 and once again in Net2.\nof partial differential equations with application to fast inverse design of nanoscale heat transport. Phys. Rev.\nResearch 4, 023210 (2022).\n(2) On the other hand, in this work, treating 𝑆 as the output, the feature of 𝐺 is not fully extracted\n𝐸\nin Net1. As shown in Fig. R19 a-b, geometry G1 and G2 share the same value of 𝑆 at about\n𝐸\n9.9 GHz, but have totally different 𝑆 . In this situation, with the same input of 𝑆 and 𝑓, 𝑆 of\n𝑅 𝐸 𝑅\nG1 and G2 could not be distinguished in the output space, if the geometric input (𝐺) is not\ntreated as the input to Net2. A similar situation can be found in geometry G3 and G4 at 8 GHz,\nas shown in Fig. R19 c-d. One difference between two situations is that 𝑆 in G1 and G2 have\n𝐸\nrates of change in the same direction, however, in G3 and G4, the changing trends are different.\nFig. R19 The input-output in Net2 where outputs differ at the same frequency with the same input of 𝑆 and\n𝐸\ndifferent 𝐺. a-b The comparison between geometry G1 and G2, where 𝑆 share the same value and similar trends\n𝐸\nat about 9.9 GHz. c-d The comparison between geometry G3 and G4, where 𝑆 share the same value but have\n𝐸\ndifferent trends at about 8 GHz.\n(3) Besides, in the revised manuscript, we introduce the major net in relation to their roles. More\nspecifically,\nSince the computational cost is high when accumulating training data for predicting 𝑆 ,\n𝑀\nthe multi-fidelity DeepOnet is used to reduce the need of high-fidelity data (𝑆 ), where data\n𝑀\nof S-parameter performance without EM-thermal coupling effects (𝑆 ) are firstly accumulated\n𝐸\nto train a low-fidelity network, i.e., Net1. The output of the low-fidelity network is further\ninput into the high-fidelity network (Net2) and formed as a residual part to reduce the difficulty\nin training Net2. To sum up, given geometry parameter (𝐺) and frequency (𝑓) of a metasurface,\nthe trained multi-fidelity DeepOnet can predict 𝑆 considering multi-physics effects with\n𝑀\nreduced requirements of high-fidelity data. To validate the effects, comparisons are added in\nthe revised manuscript (Highlighted on the last two paragraphs in “EM response by multi-\nfidelity DeepONet” in Page 8).\nFurther, since obtaining thermal effects at high frequency is either computational or\nexperimental challenging, we further developed a latent-dynamic network (Net3 & Net4) to\npredict temperature performance (𝑇(𝑓)) considering EM-thermal multi-physics coupling\neffects, where Net3 encodes 𝑆 into latent space and Net4 decodes the latent space.\n𝐸\nBy using\nEuler-method based dynamic latent, the proposed network can predict\ntemperature performance at higher frequency not seen in training data, i.e., out-of-distribution\n(OOD) prediction. To validate the effects of frequency transfer, results are presented in the\nrevised manuscript (Highlighted on the last two paragraphs in “Average temperature by latent\ndynamics networks” in Page 9).\nIn the inverse method, since designing metasurface considering both EM and thermal\nfunctionalities is a challenging multi-objective problem, hybrid regularizations are introduced\nin inverse network (Net5) to boost the training process. More specifically, besides the loss\ndefined on geometry parameters, forward solvers including EM and thermal predictor and an\nanalytical model are further introduced to alleviate non-unique problem and constrain resonant\npoints, respectively. To validate multi-functional design and the effects of different constrains,\nresults are presented in the revised manuscript (Highlighted on “Inverse design by data-\nanalytical networks” in Page 10-11).\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 4) and Supplementary Information (Highlighted on Page 5-6).\nComment 5: As opposed to in-distribution generalization that most readers would be familiar with,\nOOD generalization performance could involve a different scope, thus it is important to clearly\nspecify the scope of the generalization. Throughout the manuscript, only the generalization result\nfor the band 10-14 GHz has been reported. What would be the performance at the lower frequency\nregime (i.e., for a band < 2 GHz) could be included. Also important to observe is sharing whether\nthe frequency transfer works well even beyond 14 GHz. If not, what would the deviation pattern\nlook like? Intuitively, the error trend would presumably show a strong positive correlation w.r.t.\nthe amount of frequency deviation (i.e., the larger the frequency deviation the larger the error\nwould be). This could better present the strength/weakness of the presented frequency transfer.\nResponse: Thank you for the comment. We give our explanations and revisions as follows:\n(1) In the OOD generalization with respect to operating frequency, temperature performance\n(𝑇(𝑓)) considering EM-thermal multi-physics coupling effects is predicted by the latent-\ndynamic network, in which the latent space is updated following Euler method, which is,\n𝜕𝑎\n𝑎(𝑓 ) = 𝑎(𝑓 )+ ×𝛽 ×∆𝑓,∆𝑓 = 𝑓 −𝑓 (R6)\n𝑛+1 𝑛 𝑛+1 𝑛\n𝜕𝑓\nwhere ∆𝑓 represents the step in the frequency sequence, 𝛽 is a hyperparameter, and is\n𝜕𝑎\nchanging rate of latent space.\n𝜕𝑓\nGiven an initial state of 𝑎(𝑓 ), from (R6), the latent space can\n0\nonly be updated in one direction of increasing or decreasing frequency, which indicates the\ndirection of frequency transfer.\nHowever, the simulated EM and thermal performance of a typical geometry are given in\nFig. R20, where little information can be observed within the band below 2 GHz. Thereby, the\nprediction above 2 GHz is much more important.\nFig. R20 The simulated EM and thermal performance of a typical geometry, where the frequency goes\nbelow 2 GHz.\n(2) For the frequency transfer beyond 14 GHz, as shown in Fig. R21 a, unexpected higher order\nresonances occur at the frequency above 14 GHz, which is generally not of great concern in\nexisting works. Deviations between the predicted temperature and the full wave simulation\noccur at these resonance points, as shown in Fig. R21 b. It shows the limitation of frequency\ntransfer based on the proposed method that the prediction of different physical phenomena is\nnot accurate, which could be improved by means such as transfer learning.\nHowever, considering the grating lobe effect, the application frequency of the periodic\narray is limited by its period, i.e., 𝑓 = 𝑐⁄2𝑑= 15 GHz in this case, where 𝑓 represents the\n𝑇 𝑇\nminimum frequency grid-lobe effect may occur, 𝑐 represents the speed of light, 𝑑 represents\nthe period of array. Due to grating lobe effect and high order resonances, there is no need to\ntake such wide frequency band under consideration.\nFig. R21. Limitation of frequency transfer based on the proposed method. a Higher order resonances occur\nat the frequency above 14 GHz. b The predicted temperature compared with the results by full wave\nsimulation, which deviates at position of higher order resonances.\nTo relieve your concerns, we have added the above discussion and limitations in the revised\nmanuscript (Highlighted on Page 14) and Supplementary Information (Highlighted on Page 10).\nComment 6: Regarding the architectural design and their implementation in this work, it is\nrecommended to report the details of the architectures (e.g., the number of trainable parameters),\nas well as an example result of the network training (e.g., hyperparameters; deep learning platforms\nused for the training; optimizer) in Supporting Information. The description on the training\nprocedure under “Numerical simulation and training details” is insufficient to reproduce what’s\npresented in the manuscript. While enabling reproducibility seems optional, details of the\nimplementations could be practically informative/insightful for readers who are interested in either\nimplementing the presented framework or developing their own upon the demonstration\nshowcased by this work.\nResponse: Thank you for the comment. As shown in Fig. R22 a, Net1 and Net2 are developed in\nthe form of DeepONet based on the DeepXDE library under Tensorflow 2.13, where the dimension\nof fully connected (FC) layer is 300 in both branch-net and trunk-net. There are 4 FC hidden layers\nin branch-net and trunk-net of Net1, respectively. The number of FC hidden in in branch-net and\ntrunk-net of Net2 are 3 and 2. Thereby, the number of trainable parameters in Net1 and Net2 are\nabout 543900 and 273900. Besides, both Net1 and Net2 are trained by Adam optimizer with a\nlearning rate of 0.001.\nAs shown in Fig. R22 c, the latent dynamics networks Net3 and Net4 is developed under Keras\nand Tensorflow 2.13. Both Net3 and Net4 contain 4 FC layer of 100, where the numbers of\ntrainable parameters are 31103 and 31301, respectively. And, the hyperparameter 𝛽 is set to be\n0.05. Adam optimizer with a learning rate of 0.001 is used for training Net3 and Net4.\nAs shown in Fig. R22 b, the inverse design network Net5 is also under Keras abd Tensorflow\n2.13, which has two sub-networks. In the sub-network to deal with 𝑆 , the dimension of FC hidden\n𝐸\nlayers is [120,120,60,60,20,20]. The dimension in the other sub-network is [20,20,20]. In this case,\nthe number of trainable parameters is 35584. Adam optimizer with a learning rate of 0.001 is used\nfor training Net5.\nFig. R22. Network architecture. a The structure of multi-fidelity DeepOnet. b The structure of inverse\ndesign model. c The latent dynamics networks.\nTo relieve your concerns, we have added the above information in the Supplementary\nInformation (Highlighted on Page 5-7).\nComment 7: One of the primary motivations of deep learning assisted photonics has been to\nembrace high-dimensional, (quasi-)freeform meta-atoms. Although the authors seem to have\nachieved the frequency transfer within the specified scope, the coverage of meta-atom geometry\nseems very narrow. Based on the reviewer’s understanding the pre-defined meta-atom\nrepresentation only spans the space formed by a predefined topology, fully specified by a handful\nnumber of parameters without representation learning. Spanning and exploring only this narrow\nscope are questionable for the reviewer since operator learning, including DeepOnet, is known to\nbe powerful in dealing with high-dimensional input/output fields. The authors might want to\ninclude any discussions on this point in their manuscript.\nResponse: Thank you for the comment. We give our explanations and revisions as follows:\n(1) Besides the parameterized metasurface structures, we have also tested the frequency transfer\nof proposed method on free-form meta-atoms with higher degree of freedom. More specifically,\nwe keep the framework of the prosed method unchanged but change full connected network to\nconvolutional neural network to deal with 2D patterns. As shown in Fig. R23, the extrapolated\nfrequency band is set to be 11.2 to 14 GHz with a relative error of 0.047. It is shown that the\nfrequency transfers of the proposed method still validate in higher dimension cases.\nFig. R23. Average temperature by latent dynamics networks, where the meta-atoms are free-form and high-\ndimensional and the extrapolated frequency band is 11.2 to 14 GHz.\n(2) For more complex EM-thermal coupling effects, a three-layer meta-atom is studied in this work.\nIn this case, the top and bottom layers are the same 2D pattern with high degree of freedom, as\nshown in Fig. R24, where the simulation and prediction match well.\nTo relieve your concerns, we have added the above discussion and replaced the part of\nparameterized metasurface results with the free-form results, in the revised manuscript\n(Highlighted on Page 7) and Supplementary Information (Highlighted on Page 8-9).\nFig. R24. Average temperature by latent dynamics networks with the extrapolated frequency band at 11.2\nto 14 GHz. The meta-atoms are three-layer, in which the top and bottom layers are free-form and high-\ndimensional.\nComment 8: English could be improved. The manuscript has a lot of typos that are not difficult\nto spot. To name a few:\n- Experimental experiment\n- “Take metasurface as an example” (in page 2)\n- “However, the multi-physic field coupling problems 25-27 widen the gaps between the practical\nengineering applications and numerical simulation.” (in page 2)\n- “Experimental experiment is also conducted to verify the accuracy of the framework with a\ncomparison experiment to show a significant EM-thermal effects.” (in page 3)\nMore thorough proofreading is highly recommended.\nResponse: Thank you very much for the careful check. We have corrected the mentioned typos in\nthe revised manuscript.\nComment 9: Regarding the comment:\n“Generative adversarial network (GAN) 39,40, as one of the most popular frames in inverse\ndesign these years, constructs both forward and inverse models at the same time by adversarial\nstrategy, partly alleviating the problem of non-uniqueness.”\nThe inverse model is constructed only for conditional GAN, not the vanilla GAN (Ref. 40).\nVanilla GAN involves unsupervised learning (i.e., being agnostic to the property / output), thus what\nit does at most is “pre-processing method to compress the data space”, borrowing the authors’\nexpression. It has no direct contributions to both forward modeling (supervised learning) and\ninverse design. Thus the statement, as is, could mislead the readers who are not familiar with\nGAN. It is suggested to make a clear distinction between vanilla GAN, and conditional GAN as an\nextension.\nResponse: Thank you for the comment. In the revised manuscript, we have replaced “GAN” with\n“conditional GAN (cGAN)” and distinguished the vanilla GAN (Ref. 25) with the cGAN (Ref. 26)\nin Introduction. More specially:\nDeveloped from generative adversarial network (GAN) (Ref. 25), conditional GAN (cGAN)\n(Ref. 26), as one of the most popular frames in inverse design these years, constructs both forward\nand inverse models at the same time by adversarial strategy, partly alleviating the problem of non-\nuniqueness.\nadversarial networks. Nanophotonics 8, 1255-1261 (2019).\nComment 10: “is directed by three constraints”\nAll these are loss functions (the main objective of an optimization), not constraints (requirement /\nregularization). As another statement with a similar error:\n“Furthermore, in the area of inverse design, our approach imposes hybrid constraint to regularize\nthe results, where the analytical model extracts useful a prior from the data with equations, and\nforward discriminator makes use of knowledge by build DNN, providing the paradigm of the same\ntype of problem.”\nHerein, constraints and regularization should be rigorously distinguished.\nResponse: Thank you for the comment. We give our explanations and revisions as follows:\n(1) In general, regularization in machine learning includes explicit and implicit regularization (Ref.\n64), where explicit regularization limits model complexity and improves generalization by\nmodifying the loss function. In this work, generalized “regularization” are used to describe the\nmechanisms of inverse design method, where forward solvers and an analytical model are\nintroduced into Net5 to alleviate non-unique problem and penalize solutions that do not satisfy\nanalytical equations in training process.\nrethinking generalization. Commun. ACM 64, 107-115. (2021).\n(2) In the revised manuscript, we have replaced “constraint”, which is usually achieved by\nmodifying the definition of the optimization problem, with “regularization”.\nWe hope that the explanations and revisions provided above are acceptable to you. We sincerely\nthank you once again for your valuable time and appreciate your feedback very much.\nWith thanks and best regards,\nAuthors\nReply to the Reviewer 3\nDear Reviewer,\nWe greatly appreciate the time and insight you have invested in our paper. We provide our replies\nto your comments below and highlight the changes we have made in the manuscript using\nunderlined texts.\nComment: My comments have been mostly addressed through both the additional numerical\nresults and the revisions in the manuscript. The only remaining concern that seems partially\naddressed is Comment 4, where I have questioned whether all the component networks forming\nthe proposed forward/inverse methods are truly essential. I believe this concern, in one form or\nanother, has been commonly pointed out by the other two reviewers, while no convincing answer\ngiven by the authors can be found in both the manuscript and the rebuttal letter. I agree with\nReviewer 1 to some extent in that interlocking all the pillar networks seems not leading to\n“breakthrough”. As a workaround, the authors may want to prepare some statement on the\nscope/conditions in which the proposed framework works effectively with/without trivial\nmodifications, while conventional approaches do not. Suggesting this as an optional task, I would\nrecommend this manuscript for publication.\nResponse: Thank you for the comment. We give our explanations and revisions from the following\nthree aspects:\n(1) Importance of our forward solver (The results have also been included in reply to\nReviewer 2 in previous round): To qualitatively illustrate the importance of multi-physics\ncoupling, the influences of EM-thermal coupling effects to the performance of metasurface are\nmeasured, as shown in Fig. R4 (Fig. 6e in manuscript). From the figures, after 10 mins\nirradiation in waveguide, the surface temperature in average is 342.17 K and the measured EM\nresponse is shifted to the left during irradiation. Therefore, predicting response of metasurface\nconsidering EM-thermal coupling effect is vital.\nFig. R4 Changes in EM and thermal performance under EM-thermal coupling. a Infrared image under\n10 mins irradiation of 6-7.05 GHz, 20 W microwave. b Measured S21 during irradiation.\nImportance of our inverse solver: Furthermore, the temperature rise in metasurface under\nEM-thermal coupling could be decreased by taking the thermal-performance in inverse design,\nwhich is an improvement compared to traditional method. As shown in Fig. R5 (Fig. 5c-d in\nmanuscript), the achieved performances by the results from inverse model ignoring\ntemperature information (Option 1) and considering temperature information (Option 2) are\ngiven, respectively. From figures, with similar EM response (S21) shown in Fig. R5c, the\ntemperature of Option 1 is higher than that of Option 2.\nFig. R5 The results by inverse design considering EM-thermal multi-physics compared with those without.\na Infrared image of design option 1 under 3 mins radiation. b Infrared image of design option 2 under 3\nmins radiation. c The measured EM response of option1 and option2. d The measured temperature at the\nquarter of the samples. e The measured temperature at the half of the samples.\n(2) Ablation studies of our forward solver: To quantitatively analyze the necessity of each\ncomponent of the networks, ablation studies are conducted in the proposed forward methods.\nSpecially, in the EM forward model, as shown in Fig. R6 (Fig. 2a in manuscript), Net1 and\nNet2 together comprise the multi-fidelity DeepONet, where G and f are essential as the inputs\nto trunk-net and branch-net. Tested with 106 random samples unseen in the training process,\nthe original multi-fidelity DeepONet obtains an average relative error of 0.036. Subtracting the\nlow-fidelity results (S ) in the inputs to the trunk net in Net , the trained model, i.e., single-\nE 2\nfidelity DeepONet, obtains the relative error of 0.154 in the same validation set. Also, keeping\nS at current frequency and subtracting S (f ) and S (f ), the trained model obtains a\nE E n-1 E n+1\nvalidation error of 0.049.\nFig. R6 The structure of EM forward model, which consists of Net and Net .\n1 2\nIn the thermal forward model by latent dynamics networks, as shown in Fig. R7 (Fig. 2c in\nmanuscript), several models are trained with the same dataset and the relative error of\nfrequency transfer are compared shown in Table R1. w/o Input G refers to the situation that G\nis removed from the inputs in Net4. w/o Input S refers to the situation that S is removed from\nE E\nthe inputs in Net4. w/o Input a refers to the situation that latent space a is removed from the\ninputs in Net4. w/o Euler refers to the situation that Euler method is not used for the update of\nlatent space in Net3, which means the output to Net3 is the latent space itself. From the table,\ninput G and a(f) plays important roles in the network. In the case without Input S , because\nE\ninformation of S are partly extracted by Net3, the error in the extrapolated frequency is\nE\nrelatively small but still larger than that in the original case.\nFig. R7 The structure of thermal forward model by latent dynamics networks.\nTable R1. The relative error of frequency transfers in thermal forward model\nCase Relative error Case Relative error\nOriginal 0.043 w/o Input a 0.110\nw/o Input G 0.164 w/o Euler 0.092\nw/o Input S 0.072\nE\n(3) Ablation studies of our inverse solver (The results have been included in reply to\nReviewer 2 in previous round): In the inverse design network, the roles of EM forward model,\nthermal forward model, and analytical model is illustrated by adjusting weights in the loss\nfunction. As shown in Table. R2 (Table 2 in manuscript), the relative errors in validation set\nare given under different training case, in which the weights from w to w represent the\n1 4\nweights by geometry data, analytical model, EM forward model, and temperature forward\nmodel, respectively. E1 represents the relative error on geometry. E2 represents the relative\nerror on the resonant frequencies determined by analytical model, in which E21 and E22\nrepresent the relative error on the first and second resonant frequency, respectively. E3\nrepresents the relative error on EM response predicted by EM forward model. E4 represents\nthe relative error on average temperature predicted by thermal forward model. From the table,\nthe lowest value in each evaluation index is bolded, Case5 is thought to be the most appropriate\nchoice for our inverse design problem, for it performs the best in many indexes and has very\nclose value with the best performance in other index.\nIn Case2 and Case5, E3 is 0.151 and 0.0811 respectively, which shows the role of EM\nforward model in limiting the EM response. Comparing the E2 performance in Case3 and\nCase5, the role of analytical model in limiting resonant frequency are witnessed. Besides, the\nroles of Input-T are also studied, as shown in Fig. R8 (Fig. S11 in Supplementary Information).\nFrom the figure, under the same training weights, the relative errors in the model trained\nwithout Input-T are larger than those with Input-T, especially in thermal performance. In the\ncases without Input-T, Case7’ performs better than Case5’, which shows the necessity of\nthermal forward model in inverse design without Input-T provided.\nTable R2. The relative error under cases with different weights\nCase w w w w E1 E21 E22 E3 E4\n1 2 3 4\n1 1 0 0 0 0.315 0.0624 0.0681 0.141 0.0341\n2 0.2 0.8 0 0 0.335 0.0551 0.0600 0.151 0.0327\n3 0.3 0 0.7 0 0.294 0.0547 0.0672 0.0907 0.0305\n4 0.2 0 0 0.8 0.291 0.0610 0.0680 0.140 0.0290\n5 0.15 0.5 0.35 0 0.294 0.0539 0.0597 0.0811 0.0254\n6 0.2 0.4 0 0.4 0.302 0.559 0.0587 0.112 0.0288\n7 0.1 0.4 0.2 0.3 0.310 0.0547 0.0600 0.0893 0.0257\nFig. R8. The relative errors in validation set by inverse model trained with different weights of\nregularizations, in which the numbers from left to right in parentheses represent w to w .\n1 4\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 15) and Supplementary Information (Highlighted on Page 9).\nWe hope that the explanations and revisions provided above are acceptable to you. We sincerely\nthank you once again for your valuable time and appreciate your feedback very much.\nWith thanks and best regards,\nAuthors"}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have sufficiently addressed my concerns.", "reviewer_3": "My comments have been mostly addressed through both the additional numerical results and the revisions in the manuscript. The only remaining concern that seems partially addressed is Comment 4, where I have questioned whether all the component networks forming the proposed forward/inverse methods are truly essential. I believe this concern, in one form or another, has been commonly pointed out by the other two reviewers, while no convincing answer given by the authors can be found in both the manuscript and the rebuttal letter. I agree with Reviewer 1 to some extent in that interlocking all the pillar networks seems not leading to “breakthrough”. As a workaround, the authors may want to prepare some statement on the scope/conditions in which the proposed framework works effectively with/without trivial modifications, while conventional approaches do not. Suggesting this as an optional task, I would recommend this manuscript for publication.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Reply to the Reviewer 3\nDear Reviewer,\nWe greatly appreciate the time and insight you have invested in our paper. We provide our replies\nto your comments below and highlight the changes we have made in the manuscript using\nunderlined texts.\nComment: My comments have been mostly addressed through both the additional numerical\nresults and the revisions in the manuscript. The only remaining concern that seems partially\naddressed is Comment 4, where I have questioned whether all the component networks forming\nthe proposed forward/inverse methods are truly essential. I believe this concern, in one form or\nanother, has been commonly pointed out by the other two reviewers, while no convincing answer\ngiven by the authors can be found in both the manuscript and the rebuttal letter. I agree with\nReviewer 1 to some extent in that interlocking all the pillar networks seems not leading to\n“breakthrough”. As a workaround, the authors may want to prepare some statement on the\nscope/conditions in which the proposed framework works effectively with/without trivial\nmodifications, while conventional approaches do not. Suggesting this as an optional task, I would\nrecommend this manuscript for publication.\nResponse: Thank you for the comment. We give our explanations and revisions from the following\nthree aspects:\n(1) Importance of our forward solver (The results have also been included in reply to\nReviewer 2 in previous round): To qualitatively illustrate the importance of multi-physics\ncoupling, the influences of EM-thermal coupling effects to the performance of metasurface are\nmeasured, as shown in Fig. R4 (Fig. 6e in manuscript). From the figures, after 10 mins\nirradiation in waveguide, the surface temperature in average is 342.17 K and the measured EM\nresponse is shifted to the left during irradiation. Therefore, predicting response of metasurface\nconsidering EM-thermal coupling effect is vital.\nFig. R4 Changes in EM and thermal performance under EM-thermal coupling. a Infrared image under\n10 mins irradiation of 6-7.05 GHz, 20 W microwave. b Measured S21 during irradiation.\nImportance of our inverse solver: Furthermore, the temperature rise in metasurface under\nEM-thermal coupling could be decreased by taking the thermal-performance in inverse design,\nwhich is an improvement compared to traditional method. As shown in Fig. R5 (Fig. 5c-d in\nmanuscript), the achieved performances by the results from inverse model ignoring\ntemperature information (Option 1) and considering temperature information (Option 2) are\ngiven, respectively. From figures, with similar EM response (S21) shown in Fig. R5c, the\ntemperature of Option 1 is higher than that of Option 2.\nFig. R5 The results by inverse design considering EM-thermal multi-physics compared with those without.\na Infrared image of design option 1 under 3 mins radiation. b Infrared image of design option 2 under 3\nmins radiation. c The measured EM response of option1 and option2. d The measured temperature at the\nquarter of the samples. e The measured temperature at the half of the samples.\n(2) Ablation studies of our forward solver: To quantitatively analyze the necessity of each\ncomponent of the networks, ablation studies are conducted in the proposed forward methods.\nSpecially, in the EM forward model, as shown in Fig. R6 (Fig. 2a in manuscript), Net1 and\nNet2 together comprise the multi-fidelity DeepONet, where G and f are essential as the inputs\nto trunk-net and branch-net. Tested with 106 random samples unseen in the training process,\nthe original multi-fidelity DeepONet obtains an average relative error of 0.036. Subtracting the\nlow-fidelity results (S ) in the inputs to the trunk net in Net , the trained model, i.e., single-\nE 2\nfidelity DeepONet, obtains the relative error of 0.154 in the same validation set. Also, keeping\nS at current frequency and subtracting S (f ) and S (f ), the trained model obtains a\nE E n-1 E n+1\nvalidation error of 0.049.\nFig. R6 The structure of EM forward model, which consists of Net and Net .\n1 2\nIn the thermal forward model by latent dynamics networks, as shown in Fig. R7 (Fig. 2c in\nmanuscript), several models are trained with the same dataset and the relative error of\nfrequency transfer are compared shown in Table R1. w/o Input G refers to the situation that G\nis removed from the inputs in Net4. w/o Input S refers to the situation that S is removed from\nE E\nthe inputs in Net4. w/o Input a refers to the situation that latent space a is removed from the\ninputs in Net4. w/o Euler refers to the situation that Euler method is not used for the update of\nlatent space in Net3, which means the output to Net3 is the latent space itself. From the table,\ninput G and a(f) plays important roles in the network. In the case without Input S , because\nE\ninformation of S are partly extracted by Net3, the error in the extrapolated frequency is\nE\nrelatively small but still larger than that in the original case.\nFig. R7 The structure of thermal forward model by latent dynamics networks.\nTable R1. The relative error of frequency transfers in thermal forward model\nCase Relative error Case Relative error\nOriginal 0.043 w/o Input a 0.110\nw/o Input G 0.164 w/o Euler 0.092\nw/o Input S 0.072\nE\n(3) Ablation studies of our inverse solver (The results have been included in reply to\nReviewer 2 in previous round): In the inverse design network, the roles of EM forward model,\nthermal forward model, and analytical model is illustrated by adjusting weights in the loss\nfunction. As shown in Table. R2 (Table 2 in manuscript), the relative errors in validation set\nare given under different training case, in which the weights from w to w represent the\n1 4\nweights by geometry data, analytical model, EM forward model, and temperature forward\nmodel, respectively. E1 represents the relative error on geometry. E2 represents the relative\nerror on the resonant frequencies determined by analytical model, in which E21 and E22\nrepresent the relative error on the first and second resonant frequency, respectively. E3\nrepresents the relative error on EM response predicted by EM forward model. E4 represents\nthe relative error on average temperature predicted by thermal forward model. From the table,\nthe lowest value in each evaluation index is bolded, Case5 is thought to be the most appropriate\nchoice for our inverse design problem, for it performs the best in many indexes and has very\nclose value with the best performance in other index.\nIn Case2 and Case5, E3 is 0.151 and 0.0811 respectively, which shows the role of EM\nforward model in limiting the EM response. Comparing the E2 performance in Case3 and\nCase5, the role of analytical model in limiting resonant frequency are witnessed. Besides, the\nroles of Input-T are also studied, as shown in Fig. R8 (Fig. S11 in Supplementary Information).\nFrom the figure, under the same training weights, the relative errors in the model trained\nwithout Input-T are larger than those with Input-T, especially in thermal performance. In the\ncases without Input-T, Case7’ performs better than Case5’, which shows the necessity of\nthermal forward model in inverse design without Input-T provided.\nTable R2. The relative error under cases with different weights\nCase w w w w E1 E21 E22 E3 E4\n1 2 3 4\n1 1 0 0 0 0.315 0.0624 0.0681 0.141 0.0341\n2 0.2 0.8 0 0 0.335 0.0551 0.0600 0.151 0.0327\n3 0.3 0 0.7 0 0.294 0.0547 0.0672 0.0907 0.0305\n4 0.2 0 0 0.8 0.291 0.0610 0.0680 0.140 0.0290\n5 0.15 0.5 0.35 0 0.294 0.0539 0.0597 0.0811 0.0254\n6 0.2 0.4 0 0.4 0.302 0.559 0.0587 0.112 0.0288\n7 0.1 0.4 0.2 0.3 0.310 0.0547 0.0600 0.0893 0.0257\nFig. R8. The relative errors in validation set by inverse model trained with different weights of\nregularizations, in which the numbers from left to right in parentheses represent w to w .\n1 4\nTo relieve your concerns, we have added the above discussion in the revised manuscript\n(Highlighted on Page 15) and Supplementary Information (Highlighted on Page 9).\nWe hope that the explanations and revisions provided above are acceptable to you. We sincerely\nthank you once again for your valuable time and appreciate your feedback very much.\nWith thanks and best regards,\nAuthors"}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have addressed my remaining concern. The manuscript has been improved significantly over its initial version. I would recommend the manuscript for publication in the journal.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter