This article has been accepted for publication in a future issue of this journal, but has not been

fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 1

A Generative Machine Learning-Based Approach

for Inverse Design of Multilayer Metasurfaces
Parinaz Naseri, Student Member, IEEE, and Sean V. Hum, Senior Member, IEEE

Abstract—The synthesis of a metasurface exhibiting a spe- magnetic (TM) waves from both sides. The simulation setup
cific set of desired scattering properties is a time-consuming is implemented by stipulating periodic boundary conditions
and resource-demanding process, which conventionally relies on on the perpendicular sides of the plane and two excitation
many cycles of full-wave simulations. It requires an experienced
designer to choose the number of the metallic layers, the scatterer ports on the top and bottom of the unit cell, shown in Fig. 1
shapes and dimensions, and the type and the thickness of the (a). Deep learning [5]–[10] and statistical learning approaches
separating substrates. Here, we propose a generative machine [11] based on regression can solve the forward problem,
learning (ML)-based approach to solve this one-to-many mapping as shown in Fig. 1 (b), by providing accurate predictions
and automate the inverse design of dual- and triple-layer meta- of the scattering properties of the EMMS unit cell given its
surfaces. Using this approach, it is possible to solve optimization
problems with single or more constraints by synthesizing thin physical parameters. Such prediction tools can be used as
structures composed of potentially brand-new scatterer designs, surrogate models in analysis and optimization of uniform and
in cases where the inter-layer coupling between the layers is non-uniform metasurfaces. However, the inverse problem of
non-negligible and synthesis by traditional methods becomes predicting the physical structure of the EMMS based on the
cumbersome. Various examples to provide specific magnitude desired properties, shown in Fig 1 (c), is not easy to solve.
and phase responses of x- and y-polarized scattering coefficients
across a frequency range as well as bounded responses for The reason is that this problem is a one-to-many mapping,
different metasurface applications are presented to verify the meaning that one set of desired scattering properties might be
practicality of the proposed method. potentially provided by many different EMMSs.
Index Terms—Metasurface, inverse design, generative deep By training potentially many sub-models for different re-
learning, generative model, metasurface synthesis, machine learn- gions of the solution space, it is possible to find the optimum
ing, deep learning, surrogate models. design for simple structures composed of a single class of
scatterers [12]–[13]. A “System-by-design” scheme has been
proposed that deals with the large number of degrees of
I. I NTRODUCTION
freedom in the design of reflectarrays, where the quasi-periodic

E LECTROMAGNETIC metasurfaces (EMMSs) are 2D

structures of sub-wavelength uniform or non-uniform unit
cells composed of metallic scatterers and/or dielectric sub-
assumption for each unit cell helps to reduce the number
of the optimization variables at each iteration [14]. While
the systematic approach and the idea of optimization in a
strates. They provide the ability to manipulate electromagnetic reduced solution space are critical elements of any synthesis
waves in extraordinary ways such as spectrum filtering, wave scheme, optimization iterations of each unit cell seems to
manipulation, and polarization conversion [1]. The design of be performed through only one class of scatterer shape per
an EMMS traditionally includes two steps. The first is to reflectarray design. In the general problem where the scatterer
map the tangential field transformations on the two sides on each layer of a multilayer EMMS can have various types
of the EMMS to macroscopic properties such as scattering of shapes, both approaches become complicated by potentially
parameters, susceptibilities [2], surface impedance/admittance many submodels or optimization challenges within a high-
[3], or equivalent permittivity and permeability tensors [4]. dimensional space with many “local attraction points” [14],
The second is to implement these properties using a physical respectively. More importantly, these approaches lack the
unit cell structure. While there are analytical, optimization, ability to explore the broader design space and generate new
and systematic techniques [4] to solve the former, there is plausible solutions if required.
no straightforward method to solve this step and empirical Generative machine learning (ML) techniques have the
approaches that involve ad hoc design procedures have been capability to effectively capture the underlying patterns in
mostly employed so far. This process relies on optimiza- large complex datasets and employ that knowledge to generate
tion loops of time-consuming and resource-demanding full- new structures [15]–[16]. Several machine learning (ML)
wave simulations. While experience can expedite the design techniques for image processing have been used for the synthe-
process, adequate exploration of the design space remains sis of single-layer microwave [17]–[18] and dielectric-based
challenging. optical [19]–[21] EMMSs. Among the generative models,
Scattering properties of EMMS unit cells are usually evalu- generative adversarial networks (GANs) [18]–[21] are able to
ated when they are surrounded by similar unit cells in the same produce similar examples to the ones in the training dataset
plane and excited by transverse electric (TE) or transverse by developing an appropriate generator and critic [22]. For
generating an EMMS with desired scattering coefficients using
P. Naseri and S. V. Hum are with the Edward S. Rogers Sr. Department
of Electrical and Computer Engineering, 10 King’s College Road, toronto, GANs, a training set with the same or similar properties is
Ontario, Canada, M5S3G4, email: [Link]@[Link] required. Obtaining such a dataset is more expensive and time-

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 2

Forward Problem is potentially high-dimensional where it is difficult to find the

Periodic boundary
conditions(PBC) global optimum, and yet might not include the optimum design
Physical Scattering or
Port 1 Structure of other macroscopic due to its limited domain. Therefore, an automated generative
the EMMS properties tool is required that first traverses the design space efficiently
(b) to find the global optimum solution within it and also expands
Inverse Problem the solution space if it is inadequate. Such a tool to design an
Desired scattering or
EMMS with the desired scattering properties could greatly aid
other macroscopic Optimum Physical in designing surfaces with next-generation capabilities.
Port 2 properties Structure of the EMMS
Here, we propose a machine learning-based approach to
(a) solve the inverse problem by predicting the configuration of
(c)
thin multilayer EMMS structures given the desired scattering
Fig. 1. (a) Simulation setup to obtain the scattering coefficients properties of the surface. We leverage the interlayer coupling
of a metasurface unit cell, (b) forward, and (b) inverse problem as an extra degree of freedom to match the desired properties
of EMMS design. with the minimum error. For that, the ML tool is trained by thin
EMMSs composed of known shapes of patch- and slot-based
scatterers. Using canonical structures instead of pixelated
consuming than solving the inverse design problem itself. EMMSs, we can transfer our intelligence and experience to
Another type of generative model known as variational ML models with less training data to bring to bear the fact
autoencoders (VAE) [23] consists of an encoder and a decoder. that a great deal of fundamental EMMS behavior can be
A VAE can convert high-dimensional discrete data to a low- synthesized using these structures that do not require dense
dimensional continuous space known as the latent space, description. This tool exploits the knowledge obtained from
where standard optimization techniques can be applied to the training data and explores the solution space by generating
finding optimum solutions [24], which is not easily applicable new structures according to the new behavior they offer from
in the initial form of the input. By training the VAE, the model the canonical shapes.
learns the conditional probability distribution of the input We use an electromagnetic-aware method to generate the
given the latent variables. Therefore, it is possible to generate training data. Instead of simulating tens of thousands mul-
new plausible inputs by drawing samples from the latent tilayer EMMS structures and storing their results, which is
space and feeding them to the learned distribution. The VAE the most time-consuming and resource-demanding of a ML-
has been previously employed to represent single- and dual- based method, we only simulate the hundreds of comprising
layer metasurfaces with negligible interlayer coupling using scatterers and save their general scattering matrices (GSMs)
a continuous latent space [25]. However, the optimum latent including their higher-order mode responses, that are essential
variable is found through ad hoc random selection. Hence, this factors in capturing the interlayer coupling. The results of the
approach misses an opportunity to exploit one of the important EMMSs are then obtained by cascading these GSMs [31]. The
advantages of the VAE: that an optimization technique can cascading process takes seconds instead of minutes or hours
utilize the continuous latent space to accelerate finding the op- that the simulation of an EMMS would take. Furthermore, it
timum latent variable. Moreover, since the interlayer coupling enables us to evaluate an EMMS that is a new combination of
is negligible, the method merely relies on surrogate models the known scatterers accurately without the need of simulating
to evaluate the objectives of the optimizations, which leads to the whole structure.
inaccurate predictions in case of thin EMMSs. Furthermore, This approach utilizes VAE and regression ML models
since the desired spectra of all the reflection coefficients along as the scattering property predictors to generate multilayer
with a random latent variable are fed to the decoder, it is EMMSs based on the desired scattering coefficients. We use
not clear how much the performance of the proposed tool the proposed scheme to represent the sample EMMSs in the
gets compromised if there are no constraints on some of the training database with a continuous low-dimensional latent
properties. All the aforementioned reasons leave the need for space where interpolation, exploration, and optimization can
an approach that can design thin multilayer EMMSs for a wide be easily performed. Each latent variable represents not only
variety of constraints and applications unaddressed. the physical parameters of the EMMS, but also its scattering
Due to their increased degrees of freedom, multilayer properties. Regularizing the latent space in this way greatly
EMMSs are both promising solutions and difficult to optimize impacts the success of the optimization algorithm to efficiently
for multiobjective electromagnetic problems [26]–[28]. Plac- converge to the global optimum. We evaluate the objective
ing the layers closely using thin substrates offers additional through a hybrid approach of invoking simulation-based com-
degrees of freedom in terms of the order of behavior, but at the putations when needed, and surrogate models when they are
same time, it complicates the design further. Conventionally, known to provide sufficient accuracy, to accurately guide the
designing and optimizing structures with high interlayer cou- design process.
pling is difficult. It either requires sophisticated equivalent cir- Compared to the existing literature, the contribution of this
cuit models (ECMs) (and corresponding simulations to find the work includes: resource-efficient generation of the training
values of the ECMs) [29]–[30], or brute-forcing the problem. data; consideration of the interlayer coupling in thin multilayer
In general, the solution space to design a multilayer EMMS EMMSs; optimization in the latent space rather than doing
includes different categories of scatterer shapes. This space a brute-force search to find the global optimum design; and

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 3

efficient implementation of the simulation/ cascade-in-the-loop

to fill the gaps composed of unknown new possible designs in
Neural Neural
the latent space. 𝑧𝑘
Network Network
This paper is organized as follows. The variational autoen-
Encoder Decoder
coder is described as a generative model and compared to
conventional autoencoder in Section II. The proposed approach 𝑘 ∈ 1, 𝐾
to represent multilayer metasurfaces with a low-dimensional 𝑥 𝑥
and continuous latent space, where optimization is performed (a)
to find the optimum structure, is explained in Section III. The
details about the training data of the multilayer metasurfaces 𝜇𝑘
are presented in Section IV. The potential of the proposed Sampling
Neural Neural
approach is demonstrated through various examples in Section 𝑧𝑘
Network Network
V. Full generative capability of the proposed approach is
Encoder 𝜎𝑘 Decoder
shown in Section VI. The conclusions are drawn in Section
VII. The architectures of the neural networks are detailed in
Appendix A. 𝑥 𝑘 ∈ 1, 𝐾 , 𝑧𝑘 ~𝑁(𝜇𝑘 , 𝜎𝑘 ) 𝑥

(b)
II. T HE VARIATIONAL AUTOENCODER : A G ENERATIVE
Fig. 2. (a) Conventional autoencoder and (b) variational au-
M ODEL
toencoder (VAE).
A conventional autoencoder (AE) is an encoder followed by
a decoder, shown in Fig. 2 (a). The encoder is a neural network new examples of x given new samples from the latent space
whose input and output are a data point x and a hidden rep- based on their difference to the training samples. Depending
resentation z, respectively. Each x is described by N features on initial distribution of data in x, if the dimension of the latent
and is N -dimensional; each z is described by K number of space, K, is set to be too small, important information in x
features and is K-dimensional, where K < N . The decoder is will be lost when converted to z by the encoder. Therefore,
another neural network whose input and output are the latent the decoder cannot fully reconstruct x from z. On the other
representation z and the reconstructed x, respectively. The hand, if K is chosen to be too large, it defeats the purpose of
general idea of the AE is to learn the best encoding-decoding the VAE, and the latent space becomes sparse.
scheme using an iterative optimization process. The neural net- The loss function of a VAE, denoted by LV AE , is composed
works architectures are optimized to minimize the difference of two terms for each data point x,
between the original x and the encoded-decoded value x̂, i.e.,
reconstructed x. Thus, the overall AE architecture creates a LV AE = Lrecons + LKL . (1)
“bottleneck” for data x that ensures only the main structured This loss function is minimized over all possible x in the
part of the information can go through by the encoder and be training data set to optimize the architecture of the VAE. The
reconstructed by the decoder. AEs are employed to represent first part is the “reconstruction loss”,
the data with a lower dimensional space z and remove noise
in it. Moreover, samples with similar features will be clustered Lrecons = |x − x̂|2 . (2)
together as different separated regions in the latent space z. Minimizing Lrecons helps to convert x to z and reconstruct x̂
A variational autoencoder (VAE) is similarly the combina- from z, where x and x̂ are ideally equal. Therefore, the VAE
tion of an encoder and a decoder, shown in Fig. 2 (b). However, learns to describe x with K variables instead of N variables,
the goal of the VAE is to learn a probability distribution where K < N . The second term is the “regularization term”
P (x) over a multidimensional variable x. By modelling the that tends to organize the latent space in such a way that
distribution, it is possible to draw samples from the distribution the distributions returned by the encoder are close to the
to create new plausible values of x. Fig. 2 (b) shows the normal distribution, N (0, 1). The difference between the latent
scheme of a VAE. In practice, the encoded distributions of the space distribution and the standard normal distribution can
VAE are chosen to be normal so that the encoder can be trained be expressed by the Kulback-Leibler (KL) divergence [32],
to return the mean, µk , and the variance, σk , vectors that denoted by LKL with the closed form
describe these Gaussians, where k ∈ [1, K]. The latent space, PK
LKL = k=1 KL[N (µk , σk ), N (0, 1)]
z, is created by drawing samples from these distributions. By PK (3)
learning the distribution of z described by µk and σk instead = 21 k=1 (σk + µ2k − log(σk ) − 1).
of z itself, the VAE can be used as a generative model unlike This term gives the VAE its main and distinct feature compared
the conventional autoencoders [22]. to conventional AEs, which is encoding a set of datapoints to
Once the VAE is trained, its decoder specifically can be a continuous and differentiable latent space. It means that the
used as a generative model that outputs parameters of the VAE’s latent space does not have large gaps that would exist
likelihood distribution of p(x|z). It means that during the in the latent space of the conventional autoencoder. Over such
training process, the decoder learns to reconstruct the data a space, different optimization techniques such as gradient-
x given a representation z. After the training, it can generate descent [24] can be performed efficiently.

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 4

III. P ROPOSED A PPROACH TO O PTIMIZE M ETASURFACES

In the synthesis of a metasurface, mere generation of 𝐿𝐸𝑀𝑀𝑆
different scatterer shapes is not enough and creating samples
with a specific set of properties is required. If one is to convert
the structures of a set of EMMSs to a latent space, zgeometry ,
by the scheme shown in Fig. 2, the latent variables will
only correspond to different shapes of the scatterers and not
their scattering properties [33]. Therefore, optimization over
zgeometry can be problematic since it results in minimizing an
arbitrarily shaped loss function. For example, a loss function
to optimize an EMMS,
n
1X
LEM M S = |[S]target (i) − [S]EM M S (i)|2 , (4)
n i=1
is the mean squared error difference between the scattering
properties corresponding to each zgeometry and the desired
scattering properties, and different from the one in (1). The
scattering tensor including all the scattering coefficients of
each EMMS is denoted with [S] and n is the number of the
frequency points.
Let us first consider designing a single-layer EMMS to re- Fig. 3. An arbitrarily shaped EMMS loss function for op-
alize a specific set of desired scattering coefficients, [S]target . timizing a single-layer EMMS to obtain desired scattering
The VAE in Fig. 2 is used to convert different shapes of scat- coefficients over zgeometry , obtained by the VAE in Fig. 2.
terers to a 2-dimensional continuous latent space, zgeometry ,
where K = 2 is chosen arbitrarily for this example. The Fig. 4 shows the scheme of the proposed approach. By
dimension of the latent space is arbitrarily set to be two jointly training the VAE and the predictor, we can organize
for ease of visualization in our example. This 2-dimensional the latent space such that close variables in this space not
latent space, described by zgeometry [1] and zgeometry [2], is only have similar physical shapes, but they also present similar
composed of latent variables correspond to different shapes of scattering coefficients. This continuous latent space can be
scatterers on the bottom of Fig. 3. The corresponding EMMS used for meaningful interpolation, efficient optimization, and
loss function across the latent space for different scatterer adequate exploration. Moreover, the loss function for optimiz-
shapes is shown on the top part of in Fig. 3. As it can be seen ing the EMMS, LEM M S , gets shaped in a way that the global
from this figure, the EMMS loss function can be arbitrarily optimum can be easily found.
shaped with many local minima. This makes it difficult for
optimization methods to converge to a global optimum. There
are also some gaps in the latent space not covered by the 𝜎
Encoder Latent
samples in the training set, where the optimum latent variable space
Decoder
might lie. Therefore, we first regularize a latent space accord- 𝜇
ing to both the shape and scattering properties of the EMMSs.
frequency
Furthermore, we employ a generative model to produce new
structures based on the target scattering coefficients in case
exploring the gaps in the latent space becomes advantageous.
Predictor

𝑆
A. Using the Latent Space to Represent Metasurfaces
We adapt the scheme of the regular VAE shown in Fig. 2 so Fig. 4. The proposed approach to regularize a latent space
that the variables in the latent space not only represent different where both the physical shape of the metasurfaces and their
structures of the EMMSs but also the scattering properties they scattering properties change smoothly for better optimization.
provide. To do so, we employ an extra neural network model
as a predictor alongside the VAE. The predictor solves the
forward problem shown in Fig. 1 (a) so it outputs the scattering B. Training the Proposed ML Models
coefficients, however, its input is the latent variable instead of The scheme shown in Fig. 4 is implemented using sets of
the EMMS physical parameters. Since the latent variables are fully-connected multilayer perceptrons shown in Fig. 5. The
stochastic, we use the mean value of the latent variable, µk , as loss function, LV AEpred , is used to jointly train the neural
the input of the predictor for better training. We also use the networks shown in Fig. 5, and is defined as
frequency points as an extra input to the predictor to include
the dispersive behavior of the EMMS. LV AEpred = α × Lrecons + LKL + β × Lpred . (5)

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 5

Encoder Decoder
𝑥1 𝜖~𝑁(0,1) 𝑥1
𝑎1 𝑓1 𝑓1 𝑎1
𝜎1
𝑥2 sample
𝑎2 𝑓2 𝑓2 𝑎2 𝑥2
𝜎𝐾 𝑧1
𝑧𝑘 = 𝜇𝑘 + 𝜎𝑘 ∙ 𝜖

𝑎𝑀−1 𝑓𝐿−1 𝜇1 𝑧𝑘
𝑥𝑁−1 𝑓𝐿−1 𝑎𝑀−1
𝑥𝑁−1
𝑎𝑀 𝑓𝐿 𝜇𝐾
𝑥𝑁 𝑓𝐿 𝑎𝑀
𝑥𝑁
freq

ν1 ν2 ν𝑃−1 ν𝑃

Predictor
𝜌1 𝜌2 𝜌𝑊−1 𝜌𝑊

𝑇𝐸 𝑇𝑀 | 𝑇𝐸 ) 𝑇𝑀
|S𝑚𝑛 | |S𝑚𝑛 arg(𝑆𝑚𝑛 arg(𝑆𝑚𝑛 )

Fig. 5. Implementation of the VAE plus predictor using fully-connected multilayer perceptrons (MLP), where m, n ∈ {1, 2}.

LV AEpred is the weighted sum of the reconstruction loss, the parts of each layer, respectively. Hence, the reconstruction loss
KL divergence (3), and the prediction loss. α and β are weights is calculated using (6).
and hence hyperparameters to be tuned to achieve satisfactory The scattering properties of the multilayer EMMSs are
results in the reconstruction and the prediction processes. For obtained through a process of cascading the general scattering
small values of α and β, the KL divergence is the dominant matrices (GSMs) describing its constituent scatterers [31]. The
factor in optimizing the VAE and the predictor. This leads GSMs include high-order scattering coefficients alongside the
to a latent space with a distribution close to a predefined fundamental ones. Therefore, by doing so, we are able to
supposed distribution, but the decoder and the predictor are capture the response of the EMMS including its interlayer
not able to successfully reconstruct or predict the properties, coupling. The cascading process is fast and inexpensive.
respectively. For large values of these hyperparameters, the Therefore, creating the training data in this way is both very
reconstruction and prediction losses become the dominant efficient and resource-saving since in case a new EMMS is
factors in optimization of the weights of the neural networks of composed of a scatterer, the scatterers’ GSM can be reused.
the Fig. 5. Hence, the latent space becomes sparse with large The prediction loss, Lpred , in (5) can be calculated based
gaps where the predictor and the decoder are likely to fail to on the difference between the actual scattering coefficients of
make the right decisions for the unknown areas. Moreover, the EMMSs in the training set and the predicted ones using
both hyperparamters depend on the distributions of the data different loss functions such as the mean squared error.
in both the EMMS physical structures and their properties in The weights in the neural networks shown in Fig. 5 are
the training set. These parameters are detailed in Appendix A optimized using the backpropagation algorithm [34] by com-
for the problem at hand here. puting the derivatives of the loss function (5) with respect
The reconstruction loss, Lrecons , depends on the type of to each layer’s weights. As mentioned earlier, the latent
data in the input x. If x consists of continuous variables, the space z is stochastic and formed by sampling from a normal
reconstruction loss can be expressed using (2). Otherwise, if distribution described by µk and σk . However, it is not possible
the variables composing the input have discrete values such as to differentiate through this step and update the neural net-
0 and 1 as in the case of black and white images, we denote works parameters before it in the backpropagation algorithm.
it with X. The difference between X and X̂, reconstructed Therefore, a reparameterization trick [23] can be performed
input, can be expressed as the cross-entropy loss, to sample a variable  with standard distribution and draw
samples from the intended Gaussian as in
N
X
Lrecons = − (Xj log(X̂j ) + (1 − Xj ) log(1 − X̂j )). (6)
j=1
zk = µk + σk · k , k ∈ [1, K] and  ∼ N (0, 1). (7)

Here, each EMMS is described by the images of its layers. Since there is no need for the backpropagation algorithm
Therefore, x consists of 0 and 1 for non-metallic and metallic to pass down this sampling branch, the derivatives can be

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 6

[𝑆]𝑡𝑎𝑟𝑔𝑒𝑡
scatterer shapes. In the former case, we use a fast inexpensive
or process to cascade the stored GSMs of the known scatterers
[𝑆]𝑚𝑎𝑠𝑘𝑠 and obtain the scattering properties of the EMMS. In case the
∆[𝑆] Yes Optimum new EMMS is composed of a new scatterer shape, the image of
< ε? latent Decoder
Latent
Predictor variable the scatterer is meshed using the Rao-Wilton-Glisson (RWG)
or [𝑆]𝑙𝑣
Variable
Simulator No basis functions and fed to an in-house spectral-domain periodic
Optimum EMMS method of moments (MoM)-based simulation tool employing
those basis functions. The GSM of this scatterer is calculated
update Optimization
and saved. This process is done for all the layers of the EMMS.
The GSMs of these scatterers are similarly cascaded to obtain
[S]lv . This method eliminates the need for further simulations
in order to confirm the results. Once [S]lv is obtained using
Fig. 6. Optimization in the latent space to obtain the optimum any of the three aforementioned ways, |∆[S]| can be computed
latent variable and converting it to the optimum EMMS. analogous to (4) as
n
1X
computed as usual. The neural networks are implemented and |∆[S]| = |[S]target (i) − [S]lv (i)|2 . (8)
n i=1
trained using the TensorFlow backend Keras in Python.
It should be noted that using the fast cascading process is a
small price to pay to keep the training data small. Otherwise,
C. Optimum Metasurface
to make the predictor a reliable simulator, one needs a very
Once the VAE and the predictor in Fig. 4 are trained, large amount of training data preparing the predictor for any
the latent variables will be compact representations of the change in the shape and dimension of the EMMS scatterers.
metasurfaces, representing both their physical structures and Moreover, since scatterers with new shapes can have signifi-
scattering coefficients. Therefore, based on the target scattering cantly different scattering properties, their full-wave evaluation
coefficients, [S]target , any standard optimization algorithm is necessary to avoid proposing a wrong optimum design.
can be employed to find the optimum latent variable. Then, If |∆[S]| is less than some specified criterion, e.g. , the
the optimum latent variable is input to the decoder and the latent variable under test is considered the optimum solution
physical implementation of the optimum EMMS is obtained. in the latent space. This optimum latent variable is then
These steps are illustrated in Fig. 6. fed to the decoder so that the optimum EMMS structure is
To obtain the scattering properties corresponding to each obtained. To find the global optimum latent variable in the
latent variable, [S]lv , it matters whether the latent variable is latent space, we use the particle swarm optimization (PSO)
the compact representation of a known EMMS, i.e. it exists [35]. The optimization algorithm is implemented in Python
in the training data, or not. If the latent variable is not close using PySwarms toolkit [36].
to any of the latent variables of the training EMMSs, that
means it represents a new EMMS. For such cases, one cannot IV. D UAL -L AYER AND S YMMETRIC T HREE -L AYER
rely on the predictor to obtain the accurate results for the EMMS S AMPLES
scattering coefficients of the EMMS. This is mainly because We apply the proposed approach independently to dual-layer
for thin EMMSs where inter-layer coupling is significant, and three-layer EMMSs. In case of three-layer EMMSs, we
slight changes in the scatterers’ shapes and dimensions can consider the specific case where the EMMSs have identical
result in significantly different results. The ML models are top and bottom layers stacked with Rogers 5880 dielectric
not capable capturing such relations unless they have access slabs (r = 2.2, tanδ = 0.0009). Among the scattering
to a significant amount of training data, which is not efficient properties, without loss of generality, we train the VAE and the
and desirable for EM applications. To address this issue, we predictor on the labeled data of the TE- and TM- transmission
use a three-case method to obtain [S]lv based on the similarity coefficients of the EMMSs. Therefore, EMMSs are optimized
between the EMMS decoded from the latent variable, Xlv , and to provide desired transmission coefficients. However, one can
the EMMSs in the training data, Xt . train the models using the other scattering coefficients (such
If the latent variable gets decoded to an EMMS in the as reflection coefficients), or even entire sets or subsets of
training data or a close structure in that set, the ML predictor is scattering parameters, and optimize EMMSs based on them.
used to output [S]lv . If the minimum mean squared difference We prepare two sets of training data: one for the dual-layer
between the decoded EMMS and the EMMSs in the training EMMSs and one for the three-layer EMMSs.
data is less than some threshold, min(||Xlv − Xt ||2 ) < γ, the The primitives shown in Fig. 7 with different indicated
predictor is used. γ is specified in Section IV and is based dimensions in Table I are used for creating EMMSs for
on the primitives used in the training data. A primitive is a training. Based on the range of each parameter, the number
scatterer shape corresponding to one of the canonical ones we of shapes for each shape category is also listed in Table I.
specify. However, if the latent variable gets decoded to a new We use different shapes of asymmetric scatterers along x- and
EMMS, a simulation-based approach is used. This new EMMS y-directions such as a Jerusalem cross (Fig. 7 (a)), rectangular
might be a new combination of the known scatterer shapes, patch (Fig. 7 (b)), complementary Jerusalem cross (Fig. 7 (e)),
explained in Section IV, or combination of new generated and complementary rectangular patch (Fig. 7 (f)). In addition,

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 7

some symmetric shapes such as circular slot (Fig. 7 (c)) and thick dielectric slab are used for training. This step is faster
complete ring (Fig. 7 (d)) are used as well. These shapes and less expensive than the simulation of a dual-layer EMMS.
are chosen based on experience to provide a wide variety of Training samples of the symmetric three-layer EMMSs are
scaterring properties. generated in a similar way. About 10000 three-layer samples
Shapes in Fig. 7 (a)-(d) are used for the resonator at with dielectric thickness of 0.787 mm and 6500 samples with
the air-dielectric interface. Often alternating inductive and 1.575 mm-dielectric are created. The total time to simulate
capacitive behavior is required in EMMSs with odd number the scatterers in Fig. 7 and generate the 34000 cascaded
of layers. Therefore, we extend the shape of this scatterer to structures is 4586.7 core-hours on a Xeon(R) 4110 @2.1
include complementary Jerusalem cross and complementary GHz. The training data generation was parallelized using the
rectangular patch in Fig. 7 (e)-(f). The scatterers in dual- mutliprocessing package from Python that resulted in about
and three-layer EMMSs are simulated from 15 GHz to 31 340 hours in real time, which was spent only once for all the
GHz with unit cell period of 5.3 mm. The periodic boundary inverse problems.
conditions on x- and y-sides are stipulated. The higher-order For each dual-layer and symmetric three-layer EMMS sam-
mode scattering coefficients are the key to capture the inter- ple, the 52 × 52 image-based representations of the two
layer coupling between the resonators on different layers. scatterers of the EMMS are constructed based on Fig. 8 (a) and
Therefore, the excitation is set to be the fundamental and five fed to the encoder of the VAE. Furthermore, the decoder of
higher-order modes of x- and y- directed incident waves. This the VAE produces outputs with the size (1, 5408), whose first
number of higher-order modes provides sufficiently accurate (1, 2704) elements and second (1, 2704) elements represent
results to capture the interlayer coupling for the specified the reconstructed first and second scatterers of the EMMS,
frequency range and unit cell period. Each scatterer is trans- shown in Fig. 8 (b). In case the latent variable corresponds to
lated to meshes using RWG basis functions and fed to our a generated EMMS composed of scatterers with new shape,
in-house MoM-based simulation tool. The general scattering the 52 × 52 images of the new scatterers are post-processed to
matrix (GSM) [31] including the higher-order modes of each translate them to meshes for use with RWG basis functions,
simulation is calculated and saved. shown in Fig. 8 (b). The mesh descriptions of the scatterers can
be used to evaluate their GSM using the MoM-based solver.
𝑦 𝑦 𝑦
𝑤𝐽𝐶,𝑗𝑥 The resolution of the generated EMMS scatterers is limited by
𝑙𝑃,𝑥
𝑙𝐽𝐶,𝑗𝑥 𝑙𝐽𝐶,𝑥 the size of each pixel which is about 0.1mm×0.1mm while the
𝑥 𝑟𝐶𝑆 resolution of the known primitives shown in Fig. 7 is much
𝑥 𝑥
𝑤𝐽𝐶,𝑥 𝑙𝑃,𝑦
higher owing to the fact that pre-defined meshes are used for
them.
We train two separate sets of VAE plus predictor for dual-
(a) (b) (c)
𝑦 𝑦 𝑦
and three-layer metasurfaces and obtain distinct latent spaces
𝑤𝐶𝑅 𝑤𝐽𝐶,𝑗𝑦 for them. This is because even the same certain shapes for
𝑤𝐽𝐶,𝑦 𝑙𝑃,𝑥
𝑟𝐶𝑅 the scatterers of the dual-layer and symmetric three-layer
𝑥 𝑥 𝑥 EMMSs results in different scattering parameters. Therefore,
𝑙𝐽𝐶,𝑦 𝑙𝑃,𝑦 it is quite challenging to differentiate the two categories with
𝑙𝐽𝐶,𝑗𝑦
one latent space. Moreover, since the standard thickness of the
(d) (e) (f) dielectric is a discrete value and can change the transmission
Metal Dielectric (𝜀𝑟 = 2.2) coefficients drastically, we train different models per thickness
as well. Therefore, four different latent spaces to represent
Fig. 7. Primitives used for training: (a) Jerusalem cross (JC), dual-layer EMMSs with thicknesses of 0.787 mm and 1.575
(b) rectangular patch (RP), (c) circular slot (CS), (d) complete mm and three-layer EMMSs with total thickness of 1.574
ring (CR), (e) complementary Jerusalem cross (compJC), and mm and 3.015 mm are obtained. Information about the neural
(f) complementary rectangular patch (compRP). network architectures and their training process is detailed in
Appendix A. Based on the primitives shown in Fig. 7 and their
TABLE I. Dimensions of the Primitives in Fig. 7.
dimensions in Table I, the threshold for using the predictor,
Shape Type Parameter Value (mm) Num. of Shapes
γ, is set to be 0.03.
lJC,x/y [2.2 : 0.2 : 4.0]
lJC,jx/y lJC,x/y − 1.6 mm
To optimize both the scatterers’ shapes and the EMMS
JC & compJC 100
wJC,x/y 0.4 thickness, PSO is performed in each latent space separately
wJC,jx/jy 0.45 and the best design is chosen by selecting the best global
RP & comRP lP,x/y [2.0 : 0.2 : 5.0] 256 optimum. The swarm of the PSO includes 40 particles and it
CS rCS [1.0 : 0.1 : 2.6] 17
rCR [1.4 : 0.2 : 2.6]
is run for 50 iterations for all the examples in this manuscript.
CR 49 The number of particles was selected to make sure more than
wCR [0.1 : 0.2 : 1.3]
one particle was initially positioned in each shape category
Dual-layer EMMS training samples are generated by ran- of the training set. Moreover, the number of iterations was
domly selecting two GSMs and cascading them with different chosen to be large enough to guarantee convergence with
dielectric thicknesses [31]. About 10500 samples with 0.787 the possibility of early stopping in case no improvement was
mm-thick dielectric slab and 7000 samples with 1.575 mm- achieved after 10 iterations. All choices were made based

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 8

52 pixels
Scatterer 1:
1
𝑋𝑗,(1,52×52) 52 pixels
Scatterer 1:
1
𝑋𝑗,(1,52×52)
52 pixels
𝑋𝑗 = 𝑋𝑗1 , 𝑋𝑗2
Scatterer 2: 52 pixels (1,5408)
2
𝑋𝑗,(1,52×52)
or Scatterer 2: Input to the Encoder
2
𝑋𝑗,(1,52×52)
Scatterer 1
jth-dual-layer
jth-symmetric EMMS
three-layer EMMS
(a)
52 pixels
Scatterer 1:
1
𝑋෠𝑖,(1,2704)
𝑋෠𝑖 = 𝑋෠𝑖1 , 𝑋෠𝑖2 52 pixels
Input to the
(1,5408)
Scatterer 2:
MoM-based
Output of the Decoder 2 solver
𝑋෠𝑖,(1,2704)
of dual-layer EMMS
ith- reconstructed Conversion to RWG
dual-layer EMMS mesh files

(b)
Fig. 8. (a) The EMMS representation, i.e. Xj , to the encoder of the VAE, and (b) converting the output of the decoder to the
images of an EMMS scatterers, and then RWG basis functions to be input to the solver.

on the observation of the swarm and the achieved results smaller number of samples and iterations. In Fig. 9, it can be
for different design examples. The positions of each particle, seen that higher values of K, e.g. 12, also lead to inaccurate
denoted by lvn where lv denotes the latent variable under text predictions by the decoder and the predictor due to the sparse
and n ∈ [1, 40], is updated at the (i + 1)th iteration based on latent space.
[35]
0.06
lvn (i + 1) = lvn (i) + vn (i + 1), 𝐿𝑟𝑒𝑐𝑜𝑛𝑠
L_recons 𝐿𝑝𝑟𝑒𝑑
L_pred over test set
0.05
where
0.04
vn (i + 1) = w × vn (i) + c1 × (pn − lvn ) + c2 × (pg − lvn ).
Loss

0.03
The pn is the particle’s historically best position and pg is
the swarm’s optimum. c1 and c2 are the cognitive and social 0.02
parameters respectively. They control the particle’s behavior
given two choices: (1) to follow its personal best or (2) follow 0.01
the swarm’s global best position. Overall, this dictates if the
swarm is explorative or exploitative in nature. In addition, a 0
parameter w controls the inertia of the swarm’s movement. K=4 K=8 K = 12
Setting w = 0.8, c1 = 0.3, and c2 = 0.5 worked fine for the Fig. 9. Dependence of the reconstruction and prediction loss of
problem at our hand. One might need to tune these parameters a test set of the dual-layer EMMSs with 0.787 mm thickness.
to obtain convergence toward the global optimum for their For the purpose of visualization, this 8-dimensional latent
problems as these parameters greatly impact the performance space representing the dual-layer training set for substrate
of the PSO [35]. thickness of 1.575 mm is converted into 2-dimensions, using
The dimension K of the latent space is tuned to be 8 for t-distributed stochastic neighbor embedding (t-SNE) [37] and
minimum reconstruction and prediction losses, based on Fig. is shown in Fig. 10. From Fig. 10 (a), it can be seen that
9. As mentioned earlier, small values of K, e.g. 4, result in the latent variables corresponding to EMMSs composed of the
the loss of important information in the input, which leads to same type of scatterers, shown in Fig. 7, are clustered together.
higher reconstruction losses by the decoder. On the other hand, Moreover, we can examine this latent space in relation to the
it is encouraged to keep K as low as possible because it helps scattering properties. While the amplitude and phase of the
to 1) reduce the noise or insignificant features in the input, TE- and TM-transmission coefficients collectively impact the
i.e. the features with negligible range for the whole training latent variables, let us look at the effect of the amplitude of the
data, and 2) converge to the optimum latent variable with a TE-transmission coefficients as an example. We categorize this

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 9

property based on its number of nulls and their frequencies to bands, stopbands, and/or phase difference between the re-
5 categories: no null, one null in the range [15 − 20] GHz, sponses of the EMMS for the two orthogonal polarizations.
one null in the range [20 − 25] GHz, one null in the range The proposed method can be used for this type of general
[25 − 31] GHz, and two nulls. From Fig. 10 (b), we can see synthesis as well. To do so, we use minimum, [S]min , and max-
that there is also a correspondence between the amplitude of imum, [S]max , amplitude and phase masks over frequencies to
this coefficient and the latent variables. determine how close the latent variable’s scattering properties
are to the desired ones. Therefore, |∆[S]| in (8) is adjusted as
follows:

|∆[S]|
Pn =
i=1 ([S]lv (i) − [S]min (i)) × ([S]lv (i) − [S]max (i)) (9)
+|([S]lv (i) − [S]min (i)) × ([S]lv (i) − [S]max (i))|
It is worth noticing that if [S]lv is between the indicated bounds
by [S]min and [S]max , the |∆[S]| expectedly becomes zero.
Otherwise, it gets penalized by how much the [S]lv is out
of the bounds.
|∆[S](i)|

0
 [S]min (i) ≤ [S]lv (i) ≤ [S]max (i)
= 2 × |([S]lv (i) − [S]min (i)) × ([S]lv (i) − [S]max (i))|
(a) 
elsewhere

(10)

A. Dual-Layer EMMS for Achieving a Specific Scattering

Response
This example includes asking for an EMMS that provides
specific amplitude and phase of TE- and TM-transmission
coefficients over the frequency range. We ensured that this
set of coefficients are not achievable by any of the EMMSs
used for the training of the ML models. Therefore, a new
EMMS has to be generated to achieve the desired results. For
this example, the optimization objectives are evaluated based
on |∆[S] in (8).
(b) The amplitude and phase of the target and optimum TE
Fig. 10. The latent space of the 1.575 mm dual-layer EMMSs and TM-mode transmission coefficients are shown in Fig. 12
in the training set in regard to (a) the top and bottom layer (a), where an excellent match between the two is achieved.
scatterers and (b) the number of nulls in the amplitude of the The logarithm of the difference between the transmission
TE-transmission coefficient. coefficients corresponding to the 2-dimensional t-SNE con-
verted latent variables of the training data and the desired
V. D ESIGN E XAMPLES transmission coefficients, log10 (LEM M S ), is plotted in Fig.
In this section, the potential of the proposed approach is 12 (b). This figure shows the expected smooth change of
demonstrated through various examples of dual- and three- LEM M S over the regularized latent space. Moreover, it can be
layer EMMSs. Regardless of the application, the one-time seen that the minimum value of the log10 (LEM M S ) is −1.10,
generated training sets and trained ML models for dual- and or minimum LEM M S = 7.87 × 10−2 over the training data,
three-layer EMMSs are used for the inverse design. whereas the LEM M S corresponding to the optimum design
These examples include asking for an EMMS to achieve a is only 5.71 × 10−4 . The lower value of LEM M S for the
specific scattering response over the frequency range. While optimum design compared to this value of the best design in
using the PSO to find the global optimum latent variable, the the training set attests to the effectiveness of the optimization
flowchart in Fig. 11 is used to evaluate the scattering param- in the latent space and impact of the generative ability of the
eters of the EMMS corresponding to the latent variable under proposed approach. This optimum response is provided by a
test. For these examples, we turn the scatterer generation mode newly generated dual-layer EMMS, shown in Fig. 12 (c). This
off and explore the design space including the combinations dual-layer EMMS is the stacked Jerusalem cross and complete
of the primitives shown in Fig. 7. ring separated by 1.575 mm of Rogers 5880 dielectric. The
In practice, the desired scattering properties are rarely dimensions of the Jerusalem cross are lJC,x = 3.8 mm and
described specifically over the frequency range. Instead, one lJC,y = 2.6 mm and the rest of the parameters can be found
might prefer to define the target scattering coefficients us- in Table I. The radius and width of the ring are rCR = 2.2
ing minimum/maximum tolerable values in the transmission mm and wCR = 0.1 mm, respectively.

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 10

Latent
variable : 𝑙𝑣

Decoder

Yes
Decoded EMMS: min Xlv − Xt 2
< γ? Predictor
1 2
Scatterer shape X𝑙𝑣 = [𝑋𝑙𝑣 , 𝑋𝑙𝑣 ]
Determination No

Scatterer 1 2
𝑋𝑙𝑣 = Any 𝑋𝑙𝑣 = Any
MoM-based Yes Generation
No
primitive in
No
primitive in
simulation Mode? Fig. 7? Fig. 7?
Scatterer shape 𝑆
No Yes Yes Determination
𝑙𝑣
Closest
primitive in
Fig. 7 Get 𝐺𝑆𝑀1 Get 𝐺𝑆𝑀2

cascade
operation

Fig. 11. The flowchart to determine the scattering parameters of the EMMS corresponding to the latent space under test.

and the phase and magnitude of TM-transmission coefficient

as

|∆[S](i)|

0
 |TTE,min (i)| ≤ |Tlv (i)| ≤ |TTE,max (i)|
= 2 × |(|Tlv (i)| − |TTE,min (i)|) × (|Tlv (i)| − |TTE,max (i)|)|

elsewhere

(11)
(a)
Fig. 13 (a) shows the defined minimum and maximum
masks for a transmission coefficient whose magnitude is
between 0.9 and 1.0 for 17.0 − 19.5 GHz and 26.5 − 29.5
GHz. Moreover, the stopband is in 21.5 − 24.5 GHz where
log(𝐿𝐸𝑀𝑀𝑆 )

the amplitude is less than 0.3. It can be seen that the optimum
EMMS in Fig. 13 (b), meets the constraints with a slight
mask violation in the lower band. The optimum design is a
two Jerusalem crosses separated by a 1.575 mm Rogers 5880
dielectric. The dimensions of the top and bottom Jerusalem
crosses are based on the physical parameters in Table I, where
lJC,x = lJC,y = 3.0 and lJC,x = 3.4 mm and lJC,y = 3.0
(b) (c) mm, respectively. This optimum EMMS was generated to
Fig. 12. (a) Target and optimized TE- and TM-mode transmis- match the defined masks with much less error compared to the
sion coefficients, (b) log10 (LEM M S ) over the training t-SNE samples in the training data. It is worth noting that, in fact,
converted latent variables, and (c) optimum dual-layer EMMS. meeting the specified constraints with less error is not possible
with any dual-layer EMMS composed of a single scatterer on
each layer due to the limited order of the response that can be
B. Dual-Layer Dual-Band Frequency Selective Surface (FSS)
achieved with only two layers.
A dual-layer EMMS is optimized as a TE-polarized dual-
band frequency selective surface with a stopband between the
two transmission bands. ∆[S] in (10) is only written for the C. Three-Layer Wideband Linear-to-Circular Polarizer
magnitude of the TE-transmission coefficients while there is To design a linear-to-circular polarization converter (or
no restriction on the phase of the TE-transmission coefficient simply, a polarizer), it is required that the magnitude of the

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 11

(a)

(a)

(b)

(b)
Fig. 13. Dual-layer dual-band FSS: (a) Minimum and maxi-
mum masks for magnitude of TE-mode transmission coeffi-
cient and the optimum magnitude of TE-transmission coeffi-
cient, and (b) optimum dual-layer EMMS.

TE- and TM-transmission coefficients would be ideally equal

to 1.0 over the desired frequency band, while maintaining a
±90◦ phase difference [38]. That way, if a 45◦ -tilted linear (c)
polarized incident wave excites the EMMS, the TE and TM
components of the wave are transmitted fully with ±90◦ Fig. 14. Three-layer wideband linear-to-circular polarizer:
phase difference. Therefore, the transmitted wave would be (a) Minimum and maximum amplitude and phase difference
a circularly-polarized wave. masks for TE- and TM-transmission coefficients and the opti-
mum transmission properties, (b) optimum linear-to-circular
We design a polarizer that works from 20.5 GHz to 30.5
transmission coefficients and axial ratio, and (c) optimum
GHz. Based on the constraints explained above, we define
three-layer EMMS.
two sets of masks for the magnitude and the phase difference
between the two the transmission coefficients. The amplitude
masks are defined to keep the magnitude of both transmission D. Three-Layer Single Wideband FSS
coefficient between 0.9 and 1.0 in the frequency band of
interest. The phase masks are defined to keep the phase dif- In this example, a three-layer EMMS is optimized to
ference between 80◦ and 100◦ . Three optimization objectives filter the spectrum for the TM-mode incident waves in the
are defined for the amplitude of the TE- and TM- transmission 15.0 − 31.0 frequency band such that TM-polarized waves
coefficients and the phase difference between them, analogous between 21.0 GHz and 25.5 GHz are transmitted while those
to (11). These constraints are collected in a single cost function outside this interval get reflected. Similarly, two masks, shown
with equal and unit weights for the PSO. in Fig. 15 (a), are defined to keep the TM-transmission
The masks and the optimum transmission coefficients are coefficient amplitude above 0.9 in the transmission band and
shown in Fig. 14 (a). Fig. 14 (b) shows that the linear-to- below 0.4 in the reflection bands. The optimization objec-
circular polarization conversion is successfully done by the tive on the amplitude of the TM-transmission coefficient is
optimized generated EMMS with less than 1 dB of insertion evaluated similar to (11). The amplitude of the optimum TM-
loss and an axial ratio less than 3 dB in the desired frequency transmission coefficient and the optimum design are shown in
band. The optimized three-layer EMMS composed of a rectan- Fig. 15 (a) and (b), respectively.
gular patch with lP,x = 3.2 mm and lP,y = 2.0 mm on the top The optimum EMMS is a 1.574 mm thick structure com-
and bottom layers, and a complementary Jerusalem cross with posed of two Jerusalem crosses on the top and bottom layers
the dimensions mentioned in Table I with lJC,x = 3.0 mm and a ring in the middle layer that is embedded between
and lJC,y = 2.4 mm on the middle layer. The scatterers are two 0.787 mm Rogers 5880 substrates. The dimensions of
separated by two 1.575 mm Rogers 5880 dielectric substrates. the Jerusalem cross are listed in Table I, where lJC,x = 2.4

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 12

mm and lJC,y = 3.8. The radius and width of the ring are images of the scatterers of the optimum design are converted
rCR = 2.4 mm and wCR = 0.7 mm. This EMMS was to meshes using the RWG basis functions and simulated. The
generated to match the constraints less error compared to the GSMs of both scatterers are cascaded to obtain the scattering
samples in the training set, demonstrating the effectiveness properties of the dual-layer structure. The meshed structures
of the optimization in the latent space. It is worth noting of the scatterers on the top and bottom layers of the optimum
that the third resonance at 30 GHz is important to keep the EMMS are shown in Fig. 16 (b) and (c), respectively. These
transmission amplitude below 0.4 in the higher band. This scatterers are separated by a 1.575 mm Roger 5880 dielectric
resonance is caused by the higher-order coupling between the substrate.
scatterers in the EMMS. Therefore, besides the scatterers’ The scatterer on the top layer of this optimum EMMS,
shapes and dimensions, the thin substrate and higher-order shown in Fig. 16 (b), can be considered as the interpolation
coupling had to be optimized to achieve satisfactory results, of a rectangular patch and a Jerusalem cross. The scatterer
demonstrating the utility of the proposed approach. on the bottom layer, shown in Fig. 16 (c), appears to be an
interpolation of a Jerusalem cross and the circular slot. That
is why the latent variable corresponding to the combination of
the two scatterers lies in one of the gaps of the latent space
where it is not covered by the samples of the training set.
Therefore, the continuous representation of the EMMSs in the
latent space makes it possible to easily interpolate different
structures and obtain brand-new scattering properties.

(a)

𝑦 (a) 𝑦

(b)
Fig. 15. Three-layer single wideband TM-FSS: (a) Minimum 𝑥 𝑥
and maximum masks for magnitude of TM-mode transmission
coefficient and the optimum magnitude of TM-transmission
coefficient, and (b) optimum three-layer EMMS.

(b) (c)
VI. E XPLORATION IN THE D ESIGN S PACE
Fig. 16. (a) Minimum and maximum masks for the amplitude
Using the primitives shown in Fig. 7 (a)-(d) with the di- of the TE- and TM-transmission coefficients, (b) the top, and
mensions specified in Table I, about 178000 dual-layer EMMS (c) the bottom layers of the optimum dual-layer EMMS with
samples can be generated that provide a wide variety of TE- 1.575 mm dielectric separation.
and TM-responses. However, here, we deliberately define two
arbitrary sets of amplitude masks for the desired TE- and TM-
response, shown in Fig. 16 (a), to push the algorithm to a VII. C ONCLUSION
more “generative mode” where new scatterer shapes have to A machine learning-based approach has been proposed to
be generated to meet all the indicated requirements. Therefore, solve the inverse problem of designing compact multilayer
the “scatterer generation mode” in the shown flowchart of Fig. metasurfaces based on the desired scattering properties. Using
11 is turned on. The purpose of this example is to show that a generative ML model based on a variational autoencoder,
if the desired properties are not met by stacking the known the structures and the scattering properties of the multilayer
scatterers shown in Fig. 7, new shapes of scatterers will be EMMSs in the training set are converted to a low-dimensional
generated by the proposed approach to meet the requirements. continuous latent space. In this latent space, the particle swarm
The constraints over the magnitude of the TE- and TM- optimization has been performed to find the optimum latent
transmission coefficients are weighted equally in a single cost variable and consequently the physical design of the optimum
function for the PSO. As mentioned before, in this case, the EMMS. This approach exploits the information learned from

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 13

the training set and explores the design space by interpolating with {2048, 2048, 1024, 512, 512, 256, 128, 64} neurons and
the structures and the properties of the training structures in rectified linear unit function (ReLU) activation function. The
the latent space to propose new EMMSs that meet the desired 8th-hidden layer is connected to two 8-dimensional hidden
requirements. layers as the mean and variance layers, shown in Fig. 5. The
The EMMS optimization objectives are evaluated using dimension of the latent space is optimized to be 8 for accept-
efficient method including using ML surrogate models when able reconstruction and prediction losses. The decoder has 8
they are known to provide sufficient accuracy, fast cascading hidden layers with {64, 128, 256, 512, 512, 1024, 2048, 2048}
the GSMs of the known scatterers, and full-wave simulations neurons and the ReLU activation function. The output layer of
for brand-new generated structures. Using this method, we the decoder, denoted by x̂, has 5704 neurons with the sigmoid
eliminate the need for very large amounts of training data that activation function to reconstruct the input. It is worth noting
are required for thin EMMSs where interlayer coupling and that the sigmoid activation function is used to create values
structural modifications can change the scattering properties between 0 and 1.
significantly. Moreover, since the scattering parameters of the Frequency points between 15 and 31 GHz are normalized
optimum EMMS are always reliably computed by the tool, to values between 0 and 1 for better training and used along
there is no need for further simulation to confirm the results. the 8-dimensional latent variables as inputs of the predictor.
Therefore, while providing reliable assessment of the structure, Here, for simpler implementation, we use four predictors
we expedite the evaluation process as much as possible with to output the amplitude and phase of the TE- and TM-
the ML surrogate models and cascading process, and perform transmission coefficients, named magPredictor and phasePre-
full-wave simulations only if the proposed EMMS is likely to dictor, respectively. The magPredictors have 7 hidden layers
meet the requirements. with {500, 1000, 2000, 1000, 500, 200, 100} neurons and the
We were able not only to remove brute-forcing the combi- ReLU activation function. Since the amplitude of transmis-
nations of the known scatterers by using the PSO through a sion coefficient has a value between 0 and 1 at each fre-
latent space, but also to explore the design space by generating quency, the sigmoid activation function is used for the one-
new scatterer shapes. The optimization through the continuous dimensional output layer. The phasePredictors have 8 hidden
representation of the EMMSs greatly accelerates finding the layers with {100, 200, 500, 1000, 1000, 500, 200, 100} neurons
global optimum structure compared to using the PSO or ge- and the ReLU activation function that outputs the normalized
netic algorithms in the original space of the descriptors of the transmission phase with the sigmoid activation function.
EMMSs. Different examples of optimization problems with During the training process, the mean and variance hidden
single or multiple constraints to achieve specific dispersive layers at the end of the encoder in addition to the weights
TE- and TM-responses and application-based criteria such of the MLPs, are optimized using the Adam optimizer [39]
as dual-band frequency selective surface, wideband liner-to- with the learning rate of lr = 0.0005 to minimize the loss
circular polarization, and wideband FSS by different EMMSs function in (5). Since the input is images with 0 and 1 values,
have been demonstrated. Using individual latent spaces to the reconstruction loss, Lrecons be calculated using (6). α and
represent dual-and three-layer EMMSs with different dielectric β are tuned to 10 and 20, respectively. It is worth mentioning
thicknesses, we were able to find the global optimum by that we used gradually increasing the batch size to reduce the
adjusting both the scatterers and the dielectric thicknesses. training loss further. The reconstruction and prediction losses
The proposed approach here can be extended to inverse de- over a test set by all the trained VAEs and predictors are in
sign of metasurfaces with more than two choices of scatterers order of 1 × 10−3 and 4 × 10−4 which attest to the sufficiently
including three-layer bianistropic surfaces and surfaces with large size of the training set.
more than three layers. Moreover, other macroscopic proper-
ACKNOWLEDGMENT
ties of metasufaces such as surface admittance/impedance and
susceptibilities can be defined as desired targets based on the The authors would like to thank Zhengzheng Wang for his
application at hand. help with developing techniques to accelerate our in-house
MoM code to expedite the simulations required for this project.
A PPENDIX A R EFERENCES
N EURAL N ETWORK D ETAILS [1] O. Quevedo-Teruel et al, “Roadmap on metasurfaces,” Journal of Optics,
vol. 21, no. 7, pp. 073002 (44pp), Aug. 2019.
The representation of the EMMSs and their transmission [2] K. Achouri, M. A. Salem MA, C. Caloz, “General metasurface synthesis
coefficients can be turned to a low-dimensional latent space. based on susceptibility tensors,” IEEE Trans Antennas Propag, vol. 63,
For that, we use a jointly trained VAE and predictor, shown no. 7, pp. 2977-2991, Jul. 2015.
[3] A. Epstein and G. V. Eleftheriades, “Arbitrary power-conserving field
in Fig. 5. The two 52 × 52-images of the resonators of the transformations with passive lossless omega-type bianisotropic metasur-
EMMS are flattened and concatenated together to form a faces,” IEEE Trans. Antennas Propag., vol. 64, no. 9, pp. 3880–3895,
1 × 5408-vector. This vector, denoted by x, has 0 (no metal) 2016.
[4] G. Oliveri et al., “Synthesis of multi-layer WAIM coatings for planar
or 1 (metal) values for its components and is used as the phased arrays within the system-by-design framework,” IEEE Trans
input of the encoder. All the neural networks are implemented Antennas Propag, vol. 63, no. 6, pp. 2482-2496, Jun. 2015.
with multilayer perceptrons (MLPs). In the following para- [5] D. R. Prado, J. A. López-Fernández, G. Barquero, M. Arrebola , and
F. Las-Heras, “Fast and accurate modeling of dual-polarized reflectarray
graphs, we outline their number of hidden layers, neurons unit cells using support vector machines,” IEEE Trans. Antennas Propag.,
and activation functions. The encoder has 8 hidden layers vol. 66, no. 3, pp. 1258-1270, Mar. 2018.

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 14

[6] T. Qiu et al, “Deep learning: a rapid and efficient route to automatic [31] C. Wan and J. A. Encinar, “Efficient computation of generalized scatter-
metasurface design,”Adv. Sci., no. 6, pp. 1900128 (1-12), 2019. ing matrix for analyzing multilayered periodic structures,” IEEE Trans.
[7] V. Richard, R. Loison, R. Gillard, H. Legay, and M. Romier, “Loss Antennas Propag., vol. 43, no. 11, pp. 1233–1242, Nov. 1995.
analysis of a reflectarray cell using ANNs with accurate magnitude [32] S. Kullback and R.A. Leibler, “On information and sufficiency,” Annals
prediction,” in Proc. 11th Eur. Conf. Antennas Propag. (EuCAP), Paris, of Math. Stats., vol. 22, pp. 79–86, 1951.
France, Mar. 2017, pp. 2402–2405. [33] P. Naseri and S. V. Hum, “A machine learning-based approach to
[8] D. Kampouridou and A. Feresidis, “Machine learning-driven design synthesize multilayer metasurfaces,” IEEE International Symposium on
optimization for a multi-layer metasurface antenna,” 14th Eur. Conf. Antennas and Propagation, Montreal, QC, Canada, Jul. 2020.
Antennas Propag. (EuCAP), Copenhagen, Denmark, Mar. 2020. [34] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning internal
[9] D. Caputo, A. Pirisi, M. Mussetta, A. Freni, P. Pirinoli, and R. Zich, representations by backpropagating errors,” Nature, vol. 323, no. 6088,
“Neural network characterization of microstrip patches for reflectarray pp. 533-536, 1986.
optimization,” in Proc. 3rd Eur. Conf. Antennas Propag. (EuCAP), Berlin, [35] Kennedy and R.C. Eberhart, “Particle swarm optimization,” in Proc. of
Germany, Mar. 2009, pp. 2520-2522. the IEEE International Joint Conference on Neural Networks, pp. 1942-
[10] P. Robustillo, J. Zapata, J. A. Encinar, and J. Rubio, “ANN charac- 1948, Nov. 1995.
terization of multi-layer reflectarray elements for contoured-beam space [36] L. J. Miranda, “PySwarms: a research toolkit for Particle Swarm
antennas in the Ku-band,” IEEE Trans. Antennas Propag., vol. 60, no. 7, Optimization in Python,” Journal of Open Source Software, 3(21), 433,
pp. 3205-3214, Jul. 2012. 2018.
[11] M. Salucci et al., “Efficient prediction of the EM response of reflectarray [37] L.J.P. van der Maaten and G.E. Hinton, “Visualizing high-dimensional
antenna elements by an advanced statistical learning method,” IEEE data using t-SNE,” Journal of Machine Learning Research, no. 9,
Trans. Antennas Propag., vol. 66, no. 8, pp. 3995-4007, Aug. 2018. pp.2579-2605, Nov. 2008.
[38] P. Naseri, S. A. Matos, J. R. Costa, C. A. Fernandes and N. J. G. Fonseca,
[12] G. Gosal, “The use of inverse neural networks in the fast design of
“Dual-band dual-linear-to-circular polarization converter in transmission
printed lens antennas,” [Link]. Thesis, University of Ottawa, Ottawa,
mode application to K/Ka-band satellite communications, ” IEEE Trans.
ON, Canada, 2015, [Link]
Antennas Propag., vol. 66, no. 12, pp. 7128-7137, Dec. 2018.
[13] G. Gosal, E. Aljamali, D. McNamara, and M. Yagoub, “Transmitarray [39] D. P. Kingma and J. L. Ba, “Adam: a method for stochastic opti-
antenna design using forward and inverse neural network modeling,” mization,” International Conference on Learning Representations (ICLR),
IEEE Ant. Wireless Prop. Lett., vol. 15, no. 1, pp. 1483-1486, 2016. pages 1–13, San Diego, CA, USA, May 2015.
[14] G. Oliveri et al., “System-by-design multi-scale synthesis of task-
oriented reflectarrays,” IEEE Trans. Antennas Propag., vol. 68, no. 4,
pp. 2867-2882, Apr. 2020.
[15] Juhwan Noh et al, “Inverse design of solid-state materials via a contin-
uous representation,” Matter, vol. 1, pp.1370-1384, Nov. 2019.
[16] L. Yu, W. Zhang, J. Wang, and Y. Yu, “Seqgan: Sequence generative
adversarial nets with policy gradient,” AAAI, pp. 2852–2858, 2017.
[17] Z. Liu, D. Zhu, S. P. Rodrigues, K-T. Lee, and W. Cai, “A generative
model for inverse design of metasurfaces,” Nano Lett., vol. 18, no. 10,
pp. 6570-6576, Sep. 2018.
[18] X. Shi1, T. Qiu, J. Wang, X. Zhao, and S. Qu, “Metasurface inverse de-
sign using machine learning approaches,” Journal of Physics D: Applied
Physics, vol. 53, no.27, pp. 275105 (7pp), 2020.
[19] J. Jiang, D. Sell and J, A. Fan, “High efficiency metasurface design
based on deep generative models,” Advanced Photonics Congress (IPR,
Networks, NOMA, PVLED, SPPCom), OSA 2019.
[20] S. An et al, “Multifunctional metasurface design with a generative
adversarial network,” arXiv:1908.04851, Aug. 2019.
[21] J. Jiang and J. A. Fan, “Simulator-based training of generative models
for the inverse design of metasurfaces,” Nanophotonics, vol. 75, no. 5,
pp. 1059–1069, Nov. 2019.
[22] I. J. Goodfellow, “Generative adversarial networks,” arXiv:1406.2661.
[23] D. P. Kingma and M. Welling, “Auto-encoding variational Bayes,” In
the 2nd International Conference on Learning Representations (ICLR),
Banff, AB, Canada, Apr. 2014.
[24] R. Gomez-Bombarelli et al,“Automatic chemical design using a data-
driven continuous representation of molecules,” American Chemical So-
ciety Central Science, vol. 4, no. 2, pp. 268-276, Jan. 2018.
[25] W. Ma, F. Cheng, Y. Xu, Q. Wen, and Y. Liu, “Probabilistic represen-
tation and inverse design of metamaterials based on a deep generative
model with semi-supervised learning strategy,” Advanced Materials, vol.
31, no. 35, pp. 1901111 (9 pp.), 2019.
[26] J. Wong, M. Selvanayagam, and G. V. Eleftheriades, “A thin printed
metasurface for microwave refraction,” in Proc. IEEE MTT-S Int. Microw. Parinaz Naseri (Student Member, IEEE) received
Symp., pp. 1–4, Jun. 2014. the [Link]. and [Link]. degrees in electrical engineering
[27] H. B. Wang and Y. J. Cheng, “Single-layer dual-band linear-to-circular from the University of Tehran, Iran, in 2013 and
polarizati on converter with wide axial ratio bandwidth and different from the University of Alberta, Edmonton, AB,
polarization modes,” IEEE Trans. Antennas Propag., vol. 67, no. 6, pp. Canada, in 2017, respectively.
4296-4301, Jun. 2019. In September 2016, she received a research grant
[28] C. Zhang et al, “An ultralight and thin metasurface for radar-infrared from the Instituto de Telecomunicações, Lisbon,
bi-stealth applications”, J. Phys. D: Appl. Phys., no. 50, 444002 (7p), Portugal, to work on European Space Agency funded
2017. projects till December 2017. She joined the Recon-
[29] G. Xu, G. Elefthriades, and S. V. Hum, “Generalized synthesis technique figurable Antenna Laboratory, University of Toronto,
for high-order low-profile dual-band frequency selective surfaces,” IEEE Toronto, ON, Canada, as a Researcher in 2018,
Trans. Antennas Propag., vol. 66, no. 11, pp. 6033–6042, Nov. 2018. where she is currently pursuing the Ph.D. degree.
[30] P. Naseri, J. R. Costa, S. A. Matos, C. A. Fernandes, S. V. Hum, Ms. Naseri received Stanley G. Jones Master’s Scholarship in 2014, Ontario
“Equivalent circuit modeling to design a dual-band dual linear-to-circular Trillium Scholarship in 2018, and the Donald R. Studney Electromagnetics
polarizer surface,” IEEE Trans Antennas Propag, vol. 68, no. 7, pp.5730- Graduate Award in 2020. She is the organizer and the Chair of the IEEE UofT
5735, Jan. 2020. AP-S Student Branch Chapter.

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TAP.2021.3060142, IEEE
Transactions on Antennas and Propagation
IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. , NO. , 2020 15

Sean V. Hum (S’95–M’03–SM’11) received the

[Link]., [Link]., and Ph.D. degrees in electrical and
computer engineering from the University of Cal-
gary, Calgary, AB, Canada, in 1999, 2001, and 2006,
respectively.
In 2006, he joined The Edward S. Rogers Senior
Department of Electrical and Computer Engineering,
University of Toronto (UofT), Toronto, ON, Canada,
where he is currently a Professor and the Eugene
V. Polistuk Chair in Electromagnetic Design. He is
also the Associate Chair, Graduate Studies in the
same department. He leads the Reconfigurable Antenna Laboratory, UofT,
and along with his students, conducts research in the areas of reconfig-
urable/multifunction antennas, space-fed arrays, electromagnetic surfaces, and
antennas for space applications.
Dr. Hum was a recipient of the Governor General’s Gold Medal for his work
on radio-on-fiber systems in 2001, the IEEE Antennas and Propagation Society
Student Paper Award for his work on electronically tunable reflectarrays in
2004, the ASTech Leaders of Tomorrow Award for his work in this area in
2006, the Early Researcher Award by the Government of Ontario in 2012, and
the six UofT Departmental Teaching Awards since 2007 and a Faculty Early
Career Teaching Award in 2011 for his teaching side. He was a co-recipient
of the IEEE Antennas and Propagation Society R. W. P. King Award in 2015
and 2017. He served on the Steering Committee and the Technical Program
Committee of the 2010 and 2020 IEEE AP-S International Symposium on
Antennas and Propagation. He was the TPC Co-Chair of the 2015 IEEE
AP-S International Symposium on Antennas and Propagation. He served as
an Associate Editor for the IEEE Transactions on Antennas and Propagation
from 2010 to 2017. He is currently co-chair of the IEEE AP-S Education
Committee.

0018-926X (c) 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See [Link] for more information.
Authorized licensed use limited to: University of Cape Town. Downloaded on May 19,2021 at 20:04:45 UTC from IEEE Xplore. Restrictions apply.