# Device Interoperability for Learned Image Compression with Weights and
Activations Quantization
Esin Koyuncu12, Timofey Solovyev2, Elena Alshina2, André Kaup1

1Multimedia Communications and Signal Processing
Friedrich-Alexander-Universität Erlangen-Nürnberg
Erlangen, Germany
{andre.kaup<EMAIL_ADDRESS>2Audiovisual Laboratory, Munich Research
Center
Huawei Technologies
Munich, Germany
{elena.alshina<EMAIL_ADDRESS>
###### Abstract

Learning-based image compression has improved to a level where it can
outperform traditional image codecs such as HEVC and VVC in terms of coding
performance. In addition to good compression performance, device
interoperability is essential for a compression codec to be deployed, i.e.,
encoding and decoding on different CPUs or GPUs should be error-free and with
negligible performance reduction. In this paper, we present a method to solve
the device interoperability problem of a state-of-the-art image compression
network. We implement quantization to entropy networks which output entropy
parameters. We suggest a simple method which can ensure cross-platform
encoding and decoding, and can be implemented quickly with minor performance
deviation, of 0.3% BD-rate, from floating point model results.

###### Index Terms:

learning-based image compression, device interoperability, neural network
quantization

##  I Introduction

Image compression methods can be categorized into two approaches: traditional
codecs such as HEVC [1] and the state-of-the-art VVC [2], and learning-based
codecs such as [3], [4] and [5]. These learning-based compression networks
consist of two main levels, an analysis-synthesis part, which maps an image to
a latent representation; and a hyperprior part, which learns a prior
distribution of the latent representation for entropy coding. Prior
distributions are learned with different models in previous research, such as
zero-mean Gaussian model in [3], Gaussian model with non-zero mean such as in
[5], or Gaussian mixture model as in [4]. Also, a more complex entropy
parameter learning approach is suggested in [5] to learn Gaussian parameters
jointly, by a hyperprior and context model.

For better learning of features and representations, floating point arithmetic
is used during training of these networks. During evaluation, we need priors
to be calculated correctly at encoder and decoder side to achieve error-free
decoding. Since calculation of floating point arithmetic on different devices
is deviating, priors calculated at encoder and decoder side can be slightly
different from each other. Even a small deviation causes an error in entropy
decoding, leading to visible errors in reconstructed images at the decoder as
illustrated in Fig. 1.

Figure 1: Reconstructed images of test image 6 from decoder when encoder is
run on GPU and decoder is run on CPU. On the left: floating point model (Ours-
float). On the right: model with integerized entropy (Ours-int16-entropy).

To our knowledge, [6] is the first work addressing the learned image
compression decoder failure problem on different platforms. In their work, to
achieve determinism on different platforms, Ballé et al. propose the solution
of converting networks, weights and activations, to integer, and discretizing
the entropy parameters. A look up table (LUT) is used to calculate priors. The
authors contribute two related conclusions. The first is that if integer
networks are used for the whole model, and the model is trained with integer
emulation instead of floating point numbers, the performance of the model
degrades since using integer networks limits the approximation capacity. The
second, as a remedy for the first, is the suggestion to use integer networks
only in the prior calculation. This way, it is possible to preserve rate-
distortion performance and enable deterministic calculation. They demonstrate
results for a learned image compression model consisting of a variational
autoencoder (VAE) with a scale and zero mean hyperprior [3]. This model is an
earlier example of learned image codecs and has a 3-layer hyperdecoder to
calculate the scale entropy parameters.

Later, in [7] Sun et al. present their work for quantizing weights and
activations to 8-bit fixed point arithmetic in order to reduce complexity and
memory bandwidth, and also to ensure cross-platform decodability. In this
work, weights are quantized for both the main part and the hyper part whereas
activations are only quantized for the hyper part. As a solution to large loss
resulting from quantizing the weights of the main part, they suggest a fine-
tuning-based coding performance improvement. They also propose an extension to
the LUT of scale values, which enables sharing of CDF (cumulative distribution
function) LUTs for mean values by restricting quantized latent values to a
certain range around mean values. As a baseline, they use the model in [8]
with residual connections in the main part, a Gaussian entropy model with mean
and scale estimated by a hyperdecoder.

A more recent work to solve the platform-independent decoding problem with
post training quantization of learned image compression is reported in [9]. As
a baseline, He et al. use the image compression network model with Gaussian
mixture entropy model (GMM) estimated jointly by context and hyperprior
networks [4]. They implement integer arithmetic with a 32-bit accumulator and
quantize weights and activations to 8-bit integers only in entropy estimation
subnetworks, which results in 0.35% BD-rate degradation from the floating
point implementation for the Kodak dataset. Also, they propose an idea to
speed up the search of CDF LUTs.

In this paper we suggest a solution for the device interoperability problem by
converting parameters and activations of entropy subnetworks to 16-bit
integers. Our experiments are based on a learned image codec with GMM and
context modeling which further uses conditional coding and color separation.
Our work differs from previous work in the baseline codec being more complex,
having more layers in the entropy part, and using 16-bit for quantization of
weights and activations. The learned image compression model which we do our
experiments with is shown in Section II-A and the quantization method is
explained in Section II-B. Experimental results are presented in Section III
and conclusions are given in Section IV.

##  II Proposed method

We implement quantization to entropy subnetworks of the codec which we
experiment with. Our experiments are based on the JPEG AI Call for Proposals
(CfP)111https://jpeg.org/items/20220210_cfp_jpeg_ai.html response of TEAM14
[10] which is a variational autoencoder with color-separated conditional
coding, and with context model and GMM entropy model as shown in Fig. 2. The
CfP model consists of additional online tools providing rate distortion
optimization and rate control functionality, and post filters. However, with
the aim of a clean study of the quantization effect, we turn off all tools and
filters, i.e., use purely the neural network model in Fig. 2 in our
experiments. There are four different models trained for four different
quality ranges.

Figure 2: Architecture diagram of baseline codec, C=128.

###  II-A Description of baseline learned image codec

#### Analysis-synthesis part

In the baseline codec architecture, the luma component (Y) and chroma (U,V)
components are processed separately, and the luma component is used as an
additional input for chroma encoder and decoder. In the analysis part, there
are residual nonlinear units with attention mechanism (ResAU), simplified
residual (local) attention blocks (RAB) and gain units (GU). In Fig. 3 the
ResAU module is shown in detail. In Fig. 4 the RAB module, residual attention
module with trunk and mask branches, is shown, which is adopted from [11]. GU
provides variable rate functionality and is implemented as in [12].

Figure 3: ResAU module

Figure 4: RAB module

#### Entropy model of codec

Entropy parameters estimation is based on a hyperprior and context model. The
gather subnetwork outputs the mean, scale and weight of each Gaussian
component. We use a GMM with three Gaussian components in the entropy model.
Per color component $c=\\{Y,U,V\\}$ we can write the probability model as in
(1).

$p_{\boldsymbol{{\hat{y}}_{c}}|\boldsymbol{\hat{z}_{c}}}(\boldsymbol{{{\hat{y}}_{c}}|\boldsymbol{{\hat{z}}_{c}}})\sim\sum_{k=1}^{3}{\boldsymbol{w_{c}}^{(k)}\mathcal{N}\left(\boldsymbol{\mu_{c}}^{(k)},\boldsymbol{\sigma_{c}}^{2(k)}\right)}$
(1)

We use LeakyReLU as the activation function for all layers in the entropy
part. At the very last layer of the gather subnetwork, where we output the GMM
components, the weights, mean and scale values, Softmax is used as the
activation function. To simplify integer calculations, we calculate the
Softmax function by using the first order linear approximation of exponential
terms in the numerator and the denominator.

###  II-B Quantization method

We quantize neural network (NN) parameters, input and activations of the
neural network layers for subnetworks in the entropy part, i.e, hyper decoder,
context, and gather subnetwork, which are shown as the green boxes in Fig. 2.

The main equation for quantization to a $b$-bit integer with a restricted
range is written in (2), where $s$ is the scaling factor parameter,
clamp$(x;X_{1},X_{2})$ is the function clipping the $x$ value to the range
$[X_{1},X_{2}]$ and $\lfloor\cdot\rceil$ is rounding to the nearest integer
value. The accuracy of the quantization depends on the choice of the
parameters. In [13], authors emphasize the trade-off between clipping and
rounding error such that one can increase the clipping range in order to
decrease clipping error, but including more outliers requires a larger scaling
factor, and thus results in a larger rounding error.

${x}_{q}=\mathrm{clamp}\left(\left\lfloor\frac{x}{s}\right\rceil;-(2^{b-1}-1),2^{b-1}-1\right)$
(2)

Since we convert neural network parameters and layer inputs to 16-bit
integers, calculate convolutions with integer values and use 32-bit register
to store results, there is a risk of overflow. To make sure there is no
overflow, we calculate the scaling factors of weights quantization according
to weights data and input clipping values. We first decide on the input bit
depth, the clipping range for input of each layer, and then calculate the
weights scaling factor based on the accumulator size and the bit depth of the
input register. The decision of input clipping thresholds is made empirically,
searching the best by evaluating the model on calibration data, i.e., a few
images, and choosing the best parameter combination per layer to minimize the
bit rate. Details of quantization for activations and NN parameters are
presented below in this section.

We denote the weight parameters of the NN with
$\mathbf{W}\in\mathbb{R}^{m\times K\times K\times n}$, where $m$ and $n$ are
the number of input and output channels, respectively, and $K$ is the kernel
size of the convolution layer. $\mathbf{x}\in\mathbb{R}^{m\times h\times w}$
is the input of a layer and $\mathbf{b}\in\mathbb{R}^{n}$ is the bias
parameter.

#### Activation quantization

The main equation for quantizing the input of a layer, $\mathbf{x}$, is given
in (3).

$\mathbf{x_{q}}=\mathrm{clamp}\left(\left\lfloor\frac{\mathbf{x}}{s_{\mathbf{x}}}\right\rceil;-(2^{N_{I}-1}-1),2^{N_{I}-1}-1\right)$
(3)

We use scaling factors that are a power of two, so we decide on the shift
value, $p$, per layer which is the corresponding power of two for the scaling
factor with the relation in (4).

$s_{\mathbf{x}}=2^{-p}$ (4)

We search quantization parameters by evaluating the model on calibration data.
The shift values for the scaling factors are

  * •

$p=8$ for all layers of context, all layers of hyper decoder and the first two
layers of gather,

  * •

$p=10$ for the layers 3-7 of the gather network.

The input bit depth $N_{I}$ values are

  * •

$N_{I}=9$ for all layers of context,

  * •

$N_{I}=16$ for all layers of hyper decoder and gather.

We fuse the scaling of the output of one layer and the input of the next layer
into one operation and execute the operations in (5)-(8) in order to calculate
the quantized output of convolutions, ${\mathbf{A}}_{\mathbf{q},j}$ .

$\mathbf{{x^{{}^{\prime}}}}=\mathrm{clamp}\left(\left\lfloor{\mathbf{x}}\right\rceil;-2^{N_{I}-1}-1,2^{N_{I}-1}-1\right)$
(5) $\mathbf{A^{{}^{\prime}}}=\mathbf{b_{q}}+\mathbf{W_{q}x^{\prime}}$ (6)
${\mathbf{A}}_{\mathbf{q},j}={\frac{\lfloor{\mathbf{A}_{j}^{{}^{\prime}}}\rceil}{s_{\mathbf{W},j}{s_{\mathbf{x}}^{{}^{\prime}}}}}$
(7) ${s_{\mathbf{x}}}^{{}^{\prime}(l)}=2^{-(p^{(l)}-p^{(l-1)})}$ (8)

#### Quantization of NN parameters

Quantization of weights and bias are executed by scaling and rounding to the
nearest integer as written in (9) and (10). We do not use clipping for the
weights or bias quantization since the scaling factor is calculated to be the
maximum possible for the given weight data. Then, we calculate convolutions
with the quantized parameters $\mathbf{W_{q}}$ and $\mathbf{b_{q}}$.

$\mathbf{W}_{\mathbf{q},j,:}=\left\lfloor\frac{\mathbf{W}_{j,:}}{{s_{\mathbf{W},j}}}\right\rceil$
(9)
${\mathbf{b}_{\mathbf{q},j}}=\left\lfloor\frac{{\mathbf{b}_{j}}}{{s_{\mathbf{W},j}}s_{{\mathbf{x}}}}\right\rceil$
(10)

After deciding on the input bit depth, $N_{I}$, per layer, the per-channel
scale factor for weights, $s_{\mathbf{W},j}$ for channel $j$, are to be
calculated. We use power of two scale factors; hence, there is a corresponding
shift value, $k_{j}$, for each scaling factor $s_{\mathbf{W},j}$ with the
relation in (11).

${s_{\mathbf{W},j}}=2^{-k_{j}}$ (11)

We consider the input bit depth and the accumulator size to calculate a
maximum shift value which allows accumulation without overflow. The input bit
depth determines the clipping threshold for an input of a layer, as in (12)
$\forall i$.

$-2^{N_{I}-1}-1\leq\mathbf{x}^{{}^{\prime}}_{i}\leq 2^{N_{I}-1}-1$ (12)

For each output channel $j$, the result of the convolution of quantized input
and parameters should be within the limits which are determined by the
accumulator register size, $N_{A}=32$, so that there is no overflow. This is
shown in (13).

$-2^{N_{A}-1}\leq\sum_{i=1}^{m}{\mathbf{W}_{\mathbf{q},{ij}}\mathbf{x}^{{}^{\prime}}_{i}+\mathbf{b}_{\mathbf{q},{j}}}\leq
2^{N_{A}-1}-1$ (13)

The bias term is added to each output channel only after multiplications of
all input channels, and weights are summed up, so the bias term can be
considered at the very end. Combining restrictions from (12) and (13), we can
write (14) for the case of zero bias.

$\sum_{i=1}^{m}{|\mathbf{W}_{\mathbf{q},{ij}}|\leq 2^{N_{A}-N_{I}}}$ (14)

If we replace $\mathbf{W}_{\mathbf{q}}$ with (9), we can see the relation with
the scaling factor as in (15).

$\sum_{i=1}^{m}{\left|\left\lfloor\frac{\mathbf{W}_{ij}}{{s_{\mathbf{W},j}}}\right\rceil\right|\leq
2^{N_{A}-N_{I}}}$ (15)

Equation (16) shows the maximum per-channel shift factor, $k_{j}$, to ensure
there is no overflow for the accumulator size $N_{A}$ and the input bit depth
$N_{I}$.

${k_{j}}=N_{A}-N_{I}-\left\lceil\mathrm{log}_{2}\left(\sum_{i=1}^{m}{|\mathbf{W}_{ij}|}\right)\right\rceil$
(16)

If there is a non-zero bias term, we update the per-channel shifts for weights
calculated in (16), and choose the limiting value according to (17).

${k_{j}}\xleftarrow{}{\mathrm{min}(N_{A}-1-p-\mathrm{max}(\lceil\mathrm{log}_{2}|\mathbf{b}_{j}|\rceil,0),k_{j})-1}$
(17)

###  II-C Comparison with previous methods

The idea of using integer arithmetic in entropy subnetworks to solve the
device interoperability problem, which was proposed by Ballé et al. in [6],
was later further used in [7], [9] and also is common to our work. As Ballé et
al. point out, since encoder and decoder share these parts of the network,
they should be calculated the same in the encoder and decoder processors[6].

One difference from the previous work in ours is the baseline codec,
particularly the architecture of the entropy subnetworks of the baseline
codec. In terms of codec architecture, compared to [6] and [7], there is the
GMM entropy modeling with joint autoregressive and hyper priors in our
baseline. Compared to [9] which also use joint autoregressive priors and GMM
entropy model, the number of layers in each entropy subnetworks in our
baseline model is larger. This increases the risk of accumulation of error
resulting from the quantization of network parameters and layer inputs.

We also compare the details of quantization implementation with the previous
related work. In our implementation we choose 16-bit integers for weights and
activations and an accumulator size of 32-bit whereas [6], [7] and [9]
implement the quantization of weights and activations to 8-bit integers with
32-bit accumulators, and in [6] it is additionally suggested to implement with
4-bit weight and activation quantization with 16-bit accumulators. Even though
quantizing to a lower number of bits can be a more error-inducing operation,
and hence a more difficult problem, when both operands of convolution are
16-bit, the risk of overflow of 32-bit register is higher, compared to the
case of operands both being 8-bit. Since numeric overflow is processed
differently on different platforms, it needs to be eliminated. Our method
contributes a guarantee that there is no overflow while doing a lighter
quantization to parameters and inputs.

We propose a 16-bit quantization of entropy networks for a complex image
compression network, where we choose weight scaling factors based on the input
bit depth to eliminate the overflow of accumulator, and thus, can achieve
device interoperability, while keeping the performance of floating point
inference.

##  III Experimental Results

Our model is implemented in PyTorch. For quantization we overwrite the
convolution layers to emulate quantization with rounding and clipping
operations. We aim to reduce the BD-rate [14] loss when we quantize the
entropy subnetworks of the codec, and to eliminate the decoding error which
occurs during cross-platform coding.

#### Training

We train our model before integerization and do not update parameters during
integerization, i.e., the integerization approach does not require any
training steps. For training we use the JPEG AI training set [15].

#### Evaluation

We evaluate the model before and after quantization, and show that the
performance drop is negligible and that there is no decoder error at any of
the test images. We evaluate the performance on the Kodak dataset [16] and
Tecnick dataset [17], which consists of 100 images with a size of 1200x1200.
We calculate bit rate, and distortion as YUV-
PSNR111https://gitlab.com/wg1/jpeg-ai/jpeg-ai-qaf, and also as RGB-PSNR to
follow the literature, to evaluate the performance of the models.

Before integerization, we can observe failed reconstructed images as
illustrated in Fig. 1. With our implementation and testing on different
devices as listed in Table I, we report that for all cases decoding is
successful, i.e., there is no failure in any of the reconstructed images, and
the BD-rates between encoder and decoder, which are calculated from encoder
and decoder reconstructed images, are 0.0% for all test images.

We test our model for both float and integer-emulated versions. Table II shows
the BD-rate results, and Fig. 5 shows the corresponding rate distortion curves
for Kodak dataset, with Y-PSNR as the distortion metric.

In the previous related work, Ballé et al., Sun et al. and He et al. all
report that with their method of quantizing entropy subnetworks, they can
achieve error-free decoding without significant performance degradation from
floating point implementation. Following how the success of the previous work
is evaluated, we also show that our method does not introduce a significant
loss while providing device interoperability.

TABLE I: Tests for device interoperability Encoder device | Decoder device
---|---
NVIDIA Titan RTX | NVIDIA Titan RTX
Intel Core i9-10980XE
NVIDIA V100 | NVIDIA Titan RTX
Intel Core i9-10980XE
Intel Core i9-10980XE | NVIDIA Titan RTX
Intel Core i9-10980XE
TABLE II: BD-rate of integer entropy model vs floating point model Dataset | Y-PSNR | U-PSNR | V-PSNR | RGB-PSNR
---|---|---|---|---
Average on Kodak | 0.29 | 0.32 | 0.31 | 0.29
Average on Tecnick | 0.26 | 0.25 | 0.25 | 0.25

##  IV Discussion

In this paper, we quantize the entropy part of a learned image codec with a
hyper-prior, GMM, context models with three different kernel sizes and a
seven-layer gather, and achieve deterministic encoding and decoding on
different devices with negligible loss compared to the floating point model
despite the complexity of the entropy subnetworks. We propose a method for
choosing the quantizing parameters for 16-bit weights and activations
quantization such that the coding performance remains the same and overflow is
eliminated.

In future work, we plan to experiment with stronger quantization, quantization
of weights and activations to a smaller number of bits, and concentrate on
actual device implementation issues in order to speed up the decoding and to
enable decoding on devices with low computational capacity.

Figure 5: Rate distortion performance on Kodak data set

## References

  * [1] G. J. Sullivan, J.-R. Ohm, W.-J. Han, and T. Wiegand, “Overview of the high efficiency video coding (HEVC) standard,” IEEE Trans. Circuits Syst. Video Technol., vol. 22, no. 12, pp. 1649–1668, 2012.
  * [2] B. Bross, Y.-K. Wang, Y. Ye, S. Liu, J. Chen, G. J. Sullivan, and J.-R. Ohm, “Overview of the versatile video coding (VVC) standard and its applications,” IEEE Trans. Circuits Syst. Video Technol., vol. 31, no. 10, pp. 3736–3764, 2021.
  * [3] J. Ballé, D. Minnen, S. Singh, S. J. Hwang, and N. Johnston, “Variational image compression with a scale hyperprior,” in Proc. Int. Conf. on Learning Representations, 2018.
  * [4] Z. Cheng, H. Sun, M. Takeuchi, and J. Katto, “Learned image compression with discretized Gaussian mixture likelihoods and attention modules,” in Proc. of the IEEE/CVF Conf. on Comput. Vis. Pattern Recog., pp. 7939–7948, 2020.
  * [5] D. Minnen, J. Ballé, and G. D. Toderici, “Joint autoregressive and hierarchical priors for learned image compression,” Advances in neural information processing systems, vol. 31, 2018.
  * [6] J. Ballé, N. Johnston, and D. Minnen, “Integer networks for data compression with latent-variable models,” in Proc. Int. Conf. on Learning Representations, 2018.
  * [7] H. Sun, L. Yu, and J. Katto, “Learned image compression with fixed-point arithmetic,” in 2021 Picture Coding Symposium (PCS), (Bristol, United Kingdom), pp. 1–5, June 2021.
  * [8] Z. Cheng, H. Sun, M. Takeuchi, and J. Katto, “Deep residual learning for image compression.,” in CVPR Workshops, 2019.
  * [9] D. He, Z. Yang, Y. Chen, Q. Zhang, H. Qin, and Y. Wang, “Post-training quantization for cross-platform learned image compression,” arXiv preprint arXiv:2202.07513, 2022.
  * [10] E. Alshina, A. Boev, Z. Cui, Y. Feng, G. Gaikov, T. Guo, P. Jia, A. Karabutov, A. B. Koyuncu, E. Koyuncu, M. Li, J. Mao, S. Qian, J. Sauer, Y. Shi, T. Solovyev, M. Sychev, J. Wang, D. Yu, and Y. Zhao, “Presentation of the Huawei response to the JPEG AI Call for Proposals: Device agnostic learnable image coding using primary component extraction and conditional coding,” Tech. Rep. wg1m96016-REQ, Huawei Technologies, 2022.
  * [11] Y. Zhang, K. Li, B. Zhong, and Y. Fu, “Residual non-local attention networks for image restoration,” in Int. Conf. on Learning Representations, 2019\.
  * [12] Z. Cui, J. Wang, B. Bai, T. Guo, and Y. Feng, “G-VAE: A continuously variable rate deep image compression framework,” arXiv preprint arXiv:2003.02012, 2020.
  * [13] M. Nagel, M. Fournarakis, R. A. Amjad, Y. Bondarenko, M. van Baalen, and T. Blankevoort, “A white paper on neural network quantization,” arXiv preprint arXiv:2106.08295, 2021.
  * [14] G. Bjøntegaard, “Calculation of average PSNR differences between RD-curves,” VCEG-M33, 2001.
  * [15] “JPEG AI dataset,” 2022.  [Online]. Available ”https://jpeg.org/jpegai/dataset.html”.
  * [16] E. Kodak, “Kodak lossless true color image suite (PhotoCD PCD0992),” 1993\.  [Online]. Available ”http://r0k.us/graphics/kodak/”.
  * [17] N. Asuni and A. Giachetti, “TESTIMAGES: A large data archive for display and algorithm testing,” Journal of Graphics Tools, vol. 17, no. 4, pp. 113–125, 2013.


# Tunable Perfect Anomalous Reflection Using Passive Aperiodic Gratings

Yongming Li, Xikui Ma, Xuchen Wang, Grigorii Ptitcyn, Mostafa Movahediqomi,
and Sergei A. Tretyakov This work was supported by the China Scholarship
Council under Grant 202106280229. (Corresponding author: Xuchen Wang.)Yongming
Li is with the State Key Laboratory of Electrical Insulation and Power
Equipment, School of Electrical Engineering, Xi’an Jiaotong University, Xi’an
710049, China, and with the Department of Electronics and Nanoengineering,
Aalto University, P.O. Box 15500, FI-00076, Espoo, Finland.Xikui Ma is with
the State Key Laboratory of Electrical Insulation and Power Equipment, School
of Electrical Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049,
China.Xuchen Wang, Grigorii Ptitcyn, Mostafa Movahediqomi and Sergei A.
Tretyakov are with the Department of Electronics and Nanoengineering, Aalto
University, P.O. Box 15500, FI-00076, Espoo, Finland (e-mail:
xuchen.wang@kit.edu).

###### Abstract

Realizing continuous sweeping of perfect anomalous reflection in a wide
angular range has become a technical challenge. This challenge cannot be
overcome by the conventional aperiodic reflectarrays and periodic metasurfaces
or metagratings. In this paper, we investigate means to create scanning
reflectarrays for the reflection of plane waves coming from any direction into
any other direction without any parasitic scattering. The reflection angle can
be continuously adjusted by proper tuning of reactive loads of each array
element, while the geometrical period is kept constant. We conceptually study
simple canonical two-dimensional arrays formed by impedance strips above a
perfectly reflecting plane. This setup allows fully analytical solutions,
which we exploit for understanding the physical nature of parasitic scattering
and finding means to overcome fundamental limitations of conventional
reflectarray antennas. We propose to use subwavelength-spaced arrays and
optimize current distribution in $\lambda/2$-sized supercells. As a result, we
demonstrate perfect tunable reflection to any angle. Our work provides an
effective approach to design reconfigurable intelligent surfaces with
electrically tunable reflection angles.

###### Index Terms:

Metasurface, metagrating, reflectarray, local periodic approximation,
aperiodic, anomalous reflection.

##  I Introduction

Reconfigurable intelligent surfaces (RISs) have the prospect to become one of
the key technologies in the next generation of communication systems [1, 2,
3]. One of the main functionalities of RIS is to redirect an incoming wave to
an arbitrary direction, violating the usual reflection law. Such devices are
called anomalous reflectors, and they can be realized in several different
ways. The most well-known method is based on phase-gradient arrays and
metasurfaces [4, 5, 6, 7]. In that method, the reflection phase of meta-atoms
is designed to vary linearly along the reflector surface, using the locally
period approximation (LPA). Unfortunately, the efficiency of reflectors
designed using this method drops dramatically as the deflecting angle
increases, with strong parasitic reflections especially in the specular
direction. In recent years, it has been recognized that the cause of low
reflection efficiency is the impedance mismatch between the incident and
reflected waves. Perfect anomalous reflection requires optimization of
evanescent modes on the surface [8, 9, 10]. Therefore, more advanced methods
for perfect anomalous reflection have been proposed based on non-local
metasurfaces [11, 12, 13] and metagratings [14, 15, 16, 17, 18, 19, 20].

Most approaches to the design of metasurfaces and metagratings for reflecting
plane wave coming from a certain direction to a plane wave traveling in an
arbitrary direction are based on periodic structures. In this case, the
induced currents and scattered fields can be decomposed into spatial Fourier
series (Floquet-Bloch series). Only a finite number of Floquet modes can
propagate in space and reach the far zone, and the design goal is to maximize
the amplitude of the mode that propagates in the desired direction and
minimize the amplitudes of all other propagating harmonics. The period of the
metasurfaces and metagratings is determined by the angle of incidence and
angle of reflection. Hence, it is quite difficult to electronically scan the
angle of reflection, because scanning requires a continuous change of the
array period. However, this difficulty can be overcome by using reflectarrays.
The geometrical size of the unit cells of reflectarrays is fixed, usually to
$\lambda/2$. Known reflectarrays are designed based on the conventional
phased-array principle which has been used since the 60’s [4, 6]. In these
devices, the scanning is realized by tuning reactive loads connected to each
array element, with the goal to realize a linearly varying phase profile of
the induced current along the array plane. The loads are optimized using the
local periodic approximation (LPA), numerically simulating the normal-
incidence reflection coefficient from a single unit cell in an infinite
periodical lattice. Similarly to periodic metasurfaces designed by LPA, the
efficiency of reflectarray antennas drops dramatically if the relationship
between the incidence and reflection angles deviates much from the usual
reflection law. In particular, when the incidence angle is zero, the
efficiency becomes low if the desired angle of reflection is greater than
approximately $45^{\circ}$ [6, 21, 22]. Beyond this limitation, the parasitic
propagating modes carry parts of the energy to undesired directions. As
follows from the reciprocity principle, this problem is the same as the well-
known problem of parasitic specular reflection in off-set fed reflectarrays.

While for infinitely periodic metasurfaces, solutions to this problem are
known [23, 9, 8, 24, 10], parasitic scattering from reflectarrays is not fully
understood (different reasons for this phenomenon are identified by different
authors [22, 21]), and there are no known means to ensure perfect anomalous
reflection to any desired angle using reflectarrays with $\lambda/2$-spaced
cells. Here we test a hypothesis that for aperiodic reflectarrays, the optimal
loads of usual $\lambda/2$-spaced cells are complex-valued, and seek for
solutions to design reflectarrays that can perfectly reflect an incident wave
to arbitrary directions by using purely reactive loads. To this end, we
investigate analytically a 2D model of a tunable reflectarray, formed by
impedance-loaded thin wires or strips above a ground plane. Within this model,
for any desired distribution of induced currents in the strips, we can
analytically solve for the required load impedance of each element. Similar
structures were studied as periodic planar metagratings without scanning
capabilities [25, 26, 14, 15, 27, 28, 29, 16, 17, 19, 30, 31], but we consider
scanning finite-sized arrays with aperiodic loading. We use these analytical
solutions to find a way to design purely passive (lossless) and perfectly
operating reflectarrays for an arbitrary reflection angle. In addition to the
analytical solution, the scattered electric field distribution is calculated
by commercial software (COMSOL multiphysics).

We show that for a $\lambda/2$-spaced array the direct solution for required
load impedances indeed gives complex values, resembling the case of periodical
impenetrable impedance boundaries as anomalous reflectors [8]. We consider the
effect on performance created by dropping the real part of the calculated
loads. It is seen that the performance can go beyond the fundamental limit for
local phase-gradient reflectors, but the parasitic scattering is still
significant. To overcome this limitation, we next propose the use of several
controllable strips in each supercell of $\lambda/2$ size, which allows
optimization of the subwavelength current distributions over the array
elements. In particular, we study the arrays with three loaded strips per
$\lambda/2$, and show that proper optimization of the induced current in a
supercell allows us to realize perfect anomalous reflections in any direction.

The remainder of this article is organized as follows. In Section II, we
define the problem and develop the theory. In Section III, two different
examples are presented to validate the theory by numerical simulations. In
Section IV, we show how to realize perfect anomalous reflections using purely
reactive loads. In Section V presents the concluding observations.

##  II Principle and methodology

The anomalous reflection problem is illustrated in Fig. 1. A finite number of
thin conducting strips is placed along the $y$-direction and parallel to the
$x$-direction at distance $h$ above an infinite perfect electric conductor
(PEC) ground plane. The ground plane is placed at the $z=0$ plane. The strips
are equidistantly spaced, and the distance between them (the geometrical
period) is $d$. The strips are periodically loaded by bulk reactive loads, and
it is assumed that the distance between the load insertions is much smaller
than the wavelength of the incident wave. The structure is illuminated by a
plane wave at the incidence angle $\theta_{\rm i}$. The incident wave is
assumed to be a TE-polarized plane wave with $\mathbf{E}_{\rm
inc}(y,z)=E_{0}e^{-jk_{0}(\sin\theta_{\rm i}y+\cos\theta_{\rm i}z)}\hat{x}$.
The time dependence is assumed in the form of $e^{j\omega t}$. In the limit of
an infinite number of conducting strips, the goal is to engineer the induced
currents in the strips so that the incident wave would be fully reflected to
the anomalous direction at $\theta_{\rm r}$. If we want to generate only one
plane wave along the $\theta_{\rm r}$ direction, the simplest solution is to
ensure that strips support induced currents whose phase varies linearly with
the phase gradient dictated by the reflection angle $\theta_{\rm r}$ (e.g.,
[5, 28]). Indeed, according to the phased-array antenna theory, if the
currents flowing in an infinite array of strips have uniform amplitude and a
linearly varying phase profile, a plane wave is generated [32, 33]. In
addition, it is necessary to eliminate the specular reflection of the incident
wave in the PEC ground plane, in order to achieve full-power reflection in the
desired direction. For a finite number of strips, such design ensures the
highest directivity of the beam reflected towards $\theta_{\rm r}$ (equivalent
to the pattern of a uniformly fed active phased array).

In practice, thin conducting strips are modeled as equivalent round wires with
the effective radius $r_{\rm eff}=w/4$ [32]. For a single current line
$\mathbf{J}=I\delta(y_{0},z_{0})\hat{x}$ positioned at $(y_{0},z_{0})$ and
parallel to $x$-axis, the scattered electric field at point $(y,z)$ can be
calculated as [34, 32]

$\mathbf{E}_{x}=-\frac{k_{0}\eta_{0}}{4}IH_{0}^{(2)}(k_{0}\sqrt{(y-y_{0})^{2}+(z-z_{0})^{2}})\hat{x},$
(1)

where $k_{0}=\omega\sqrt{\mu_{0}\epsilon_{0}}$ is the wavenumber in free
space, $\eta_{0}=\sqrt{\mu_{0}/\epsilon_{0}}$ represents the wave impedance in
free space, and $H_{0}^{(2)}(\cdot)$ is the zeroth-order Hankel function of
the second kind. Due to the presence of the ground plane, the electric field
produced by strips is created by the induced currents flowing in the strips
and their images. The specular reflection of the incident wave from the ground
plane produces the field

$\mathbf{E}_{\rm ref}=-E_{0}e^{-jk_{0}(\sin\theta_{\rm i}y-\cos\theta_{\rm
i}z)}\hat{x}.$ (2)

The total external field that excites the strip array is the sum of the
incident and the specularly reflected waves, i.e., $\mathbf{E}^{\rm
ext}(y,z)=-j2E_{0}\sin(k_{0}\cos\theta_{\rm i}z)e^{-jk_{0}\sin\theta_{\rm
i}y}\hat{x}$.

In order to eliminate the specular reflection and generate the desired
reflected plane wave, two sets of induced currents are required. Let us denote
one set of currents that eliminates specular reflection as
$\mathbf{J}_{\alpha}(y,z)=\sum_{m=0}^{N-1}{I_{\alpha}\over
d}e^{-jk_{0}\sin\theta_{\rm i}y_{m}}\delta(y-y_{m},z+h)\hat{x}$, and the other
set of currents that creates the desired anomalously reflected wave as
$\mathbf{J}_{\beta}(y,z)=\sum_{m=0}^{N-1}{I_{\beta}\over
d}e^{-jk_{0}\sin\theta_{\rm i}y_{m}}\delta(y-y_{m},z+h)\hat{x}$. Because both
current components generate plane waves (in the limit of an infinite array,
the complex amplitudes of the strip currents $I_{\alpha,\beta}$ are the same
for all strips.

Figure 1: (a) An array of $N$ strips over an infinite ground plane under a
plane-wave illumination at $\theta_{\rm i}$. The desired reflection angle is
$\theta_{\rm r}$. The fields radiated by the currents induced in the strips
eliminate the specular reflection of the incident wave in the ground plane and
launch the desired anomalously reflected plane wave. The distance between
strips is $d$. (b) Front view of the array. The strips are loaded by bulk
impedances inserted periodically with the period $l$. The width of the strips
is $w$. Both $l$ and $w$ are much smaller than the wavelength in free space.
The periodically loaded strips can be modeled as homogeneous impedance strips
with the impedances per unit length $Z_{{\rm L},n}$, where
$n\in\\{0,1,2,\cdots,N-1\\}$.

###  II-A Calculation of current density distribution

Let us assume that the current density flowing in strip $m$ is
$\mathbf{J}_{m}(y,z)=I_{m}\delta(y-y_{m},z+h)\hat{x}$, where
$y_{m}=md~{}(m\in\\{0,1,2,\cdots,N-1\\}$ is the position of the loaded strip.
The total current $I_{m}$ in the strip $m$ has dimension Ampere. The electric
field generated by current line $m$ at the position of current line $n$
$(m\neq n)$ is

$\mathbf{E}_{nm}=-\frac{k_{0}\eta}{4}I_{m}\left[H_{0}^{(2)}(k_{0}y_{mn})-H_{0}^{(2)}(k_{0}\sqrt{y_{mn}^{2}+4h^{2}})\right]\hat{x},$
(3)

where $y_{mn}=\left|m-n\right|d$ is the distance between the strips $m$ and
$n$. The second term of (3) is the electric field generated by the mirror
image of the strip in the PEC plane. The corresponding mutual impedance
between strips $m$ and $n$ can be calculated as
$Z_{nm}=\frac{k_{0}\eta}{4}\left[H_{0}^{(2)}(k_{0}y_{mn})-H_{0}^{(2)}(k_{0}\sqrt{y_{mn}^{2}+4h^{2}})\right]$.

According to Ohm’s law, the current in the loaded strip with the index $n$
satisfies

$Z_{{\rm L},n}I_{n}=E^{\rm ext}_{x}(y_{n},-h)-\sum_{m=0,\\\ m\neq
n}^{N-1}Z_{nm}I_{m}-Z_{n}I_{n},$ (4)

where $Z_{n}={k_{0}\eta\over 4}\left[H_{0}^{(2)}(k_{0}r_{\rm
eff})-H_{0}^{(2)}(2k_{0}h)\right]$ are the self-impedances of strips per unit
length. The left-hand side is the total electric field on the strip numbered
$n$. The first term on the right-hand side is the external electric field at
the strip numbered $n$, _i.e._ , $E_{x}^{\rm
ext}(y_{n},-h)=j2E_{0}\sin(k_{0}\cos\theta_{\rm i}h)e^{-jk_{0}\sin\theta_{\rm
i}nd}$. The second term on the right-hand side is the field created due to
mutual interactions with the other strips. The last term on the right-hand
side is caused by its self-action. The current in each loaded strip has to
satisfy (4), which can be written in a matrix form

$\overline{\overline{Z}}\cdot\vec{I}=\vec{U},$ (5)

where $\vec{I}=\left[I_{0},I_{1},\cdots,I_{N-1}\right]^{\rm T}$ is the vector
of induced currents and $\vec{U}=\left[E^{\rm ext}_{x}(y_{0},-h),E^{\rm
ext}_{x}(y_{1},-h),\cdots,E^{\rm ext}_{x}(y_{N-1},-h)\right]^{\rm T}$ is the
voltage vector. $\overline{\overline{Z}}=\overline{\overline{Z}}_{\rm
s}+\overline{\overline{Z}}_{\rm m}+\overline{\overline{Z}}_{\rm L}$ is the
impedance matrix that is composed of the self-impedance matrix (a diagonal
matrix) $\overline{\overline{Z}}_{\rm s}={\rm
diag}\left(Z_{0},Z_{1},\cdots,Z_{n},\cdots,Z_{N-1}\right)$, the load impedance
matrix $\overline{\overline{Z}}_{\rm L}={\rm diag}(Z_{{\rm L},0},Z_{{\rm
L},1},\cdots,Z_{{\rm L},N-1})$, and the mutual impedance matrix
$\overline{\overline{Z}}_{\rm m}$ (see (6)), which has the form

$\overline{\overline{Z}}_{\rm
m}=\left(\begin{array}[]{ccccc}0&Z_{0,1}&\cdots&Z_{0,N-1}\\\
Z_{1,0}&0&\cdots&Z_{1,N-1}\\\ Z_{2,0}&Z_{2,1}&\cdots&Z_{2,N-1}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
Z_{N-1,0}&Z_{N-1,1}&\cdots&0\end{array}\right).$ (6)

As discussed above, the required induced currents flowing in the conductive
strips should contain two components, one launches a plane wave to eliminate
the specular reflection from the ground plane, and the other generates the
desired plane wave. The amplitudes of these two components should be uniform
over the whole array. We denote the amplitudes of the first current component
as $I_{\alpha}$, and the second component as $I_{\beta}$. The phases of these
two components are set as linear functions along the array plane, to ensure
launching waves into the required directions:
$I_{m}=I_{\alpha}e^{-jk_{0}\sin\theta_{\rm
i}y_{m}}+I_{\beta}e^{-jk_{0}\sin\theta_{\rm r}y_{m}}$. Next, we will show how
to calculate the required induced currents amplitudes.

####  II-A1 Elimination of the specularly reflected wave

To eliminate specular reflection from the ground plane, the sum of the
surface-averaged value (the fundamental Floquet harmonic) of the scattered
field generated by the induced strip currents and the incident field reflected
from the ground should equal zero. We write this condition at $z=-h$, just
above the array: $\mathbf{E}_{\alpha}+\mathbf{E}_{\rm ref}(y,-h)=0$. The
electric field generated by the strips can be found using [32, Eq. 4.35],
which gives $-\frac{\mathbf{J}_{\alpha}\eta_{0}}{2\cos\theta_{\rm i}}$. The
total radiated field is the sum of that electric field generated by the strips
and the electric field reflected from the ground:

$\mathbf{E}_{\alpha}=-\frac{\mathbf{J}_{\alpha}\eta_{0}}{2\cos\theta_{\rm
i}}\left(1-e^{-jk_{0}\cos\theta_{\rm i}2h}\right).$ (7)

Substituting $\mathbf{E}_{\rm ref}(y,-h)$ from (2) and writing
$\mathbf{J}_{\alpha}={I_{\alpha}\over d}\hat{x}$, we find

$I_{\alpha}=j\frac{E_{0}d\cos\theta_{\rm i}}{\eta_{0}\sin(k_{0}\cos\theta_{\rm
i}h)}.$ (8)

####  II-A2 Generation of the desired plane wave

The averaged scattered electric field generated by the current component
$I_{\beta}$ reads, similarly,

$\mathbf{E}_{\beta}=-\frac{\mathbf{J}_{\beta}\eta_{0}}{2\cos\theta_{\rm
r}}(1-e^{-jk_{0}\cos\theta_{\rm r}2h}).$ (9)

From this relation we see that the distance between the plane of strips and
the ground plane should satisfy $h\neq l\frac{\lambda_{0}}{2\cos\theta_{\rm
r}}$, where $l$ is a positive natural number. Otherwise, the sum of the
primary field of these currents and the reflection from the ground would equal
zero, $\mathbf{E}_{\beta}=0$, and we could never generate the desired plane
wave.

To find the required amplitude of current $\mathbf{J}_{\beta}$, we will use
the power conservation, equating the normal components of the incident plane
wave and the anomalously reflected plane wave. First, we calculate the
tangential component of the magnetic field that corresponds to the plane wave
with the electric field (9). This gives

$\mathbf{H}_{\beta}\cdot\hat{y}=\frac{E_{\beta}}{\eta_{0}}\cos\theta_{\rm r}.$
(10)

The corresponding normal component of the Poynting vector reads
$\frac{|E_{\beta}|^{2}}{2\eta_{0}}\cos\theta_{\rm r}$. Equating that to the
normal component of the incident-wave Poynting vector $\frac{|E_{\rm
inc}|^{2}}{2\eta_{0}}\cos\theta_{\rm i}$, we find the required current
amplitude

$\left|{I_{\beta}}\right|=|J_{\beta}|d=\left|\frac{E_{0}d\sqrt{\cos\theta_{\rm
i}\cos\theta_{\rm r}}}{\eta_{0}\sin(k_{0}h\cos\theta_{\rm r})}\right|.$ (11)

Naturally, the phase of current $I_{\beta}$ does not affect the power carried
by the reflected plane wave. However, it does not mean that this phase has no
significance in reflector design and applications. Actually, if there is no
need to engineer a specific reflection phase, the phase of current $I_{\beta}$
can be used as an optimization variable. Later, in Section IV, we will see
that when the real parts of load impedances are neglected, the efficiency can
be improved by optimizing the phase of $I_{\beta}$.

###  II-B Efficiency

Once the induced current distribution is known, the electric field generated
by the strip array can be calculated as E^strips=&\- k0η04 ∑_m=0^N-1 I_m [
H_0^(2) (k_0 (y - y_m )^2+ (z+h )^2 )
\- H_0^(2)(k_0 (y - y_m )^2+ (z-h )^2 ) ] ^x. The scattered electric field is
calculated by subtracting the incident electric field from the total electric
field. For our finite-sized array, the scattered electric field
$\mathbf{E}^{\rm sca}$ is the sum of the field generated by the strips
$\mathbf{E}^{\rm strips}$ and the specularly reflected field $\mathbf{E}_{\rm
ref}$ (2), _i.e._ ,

$\mathbf{E}^{\rm sca}=\mathbf{E}^{\rm strips}+\mathbf{E}_{\rm ref}.$ (12)

Usually, the efficiency of anomalous reflection is defined for infinite
periodical reflectors as the ratio of the power densities carried by the
reflected and incident plane waves or amplitude of the tangential component of
the reflected wave and incident wave. For a finite-sized array, it is
impossible to generate an ideal plane wave that does not attenuate with the
direction of propagation. Since the amplitude of the scattered electric field
decreases along the propagation direction, we define the power efficiency of
anomalous reflection in terms of far-zone fields as

$\zeta=\frac{|E^{\rm far}(\theta_{\rm r},\vec{Z}_{\rm L})|^{2}}{|E_{\rm
refer}^{\rm far}(\theta_{\rm r})|^{2}},$ (13)

where $E^{\rm far}(\theta,\vec{Z}_{\rm L})$ is the electric field amplitude
generated by strips $E^{\rm strips}$ in the far zone along the direction
$\theta_{\rm r}$ for a specific set of strip load impedances which can be
written into a vector $\vec{Z}_{\rm L}=\left(Z_{{\rm L},0},Z_{{\rm
L},1},\cdots,Z_{{\rm L},N-1}\right)^{\rm T}$. $E_{\rm refer}^{\rm
far}(\theta_{\rm r})$ is the reference far-zone electric field along the same
desired direction created by the ideal current distribution of $I_{\alpha}$
and $I_{\beta}$ (with uniform amplitude and linearly varying phase) in Eq.
(II-B).

##  III Analytically calculated required load impedances

Now that we know what induced currents should flow in the wires in order to
anomalously reflect a given plane wave into the desired direction, we can
calculate the corresponding load impedances that need to be connected to the
strips. This can be done algebraically by solving the system of equations (5)
for the matrix of impedance loads ${\overline{\overline{Z}}}_{\rm L}$. Indeed,
we know the vector of incident fields $\vec{U}$, the vector of required
currents in the strips $\vec{I}$, and the impedance matrix of the strip array,
and the only unknowns are the load impedances. As a particular example, we
assume an incident wave with the frequency $f=10~{}\text{GHz}$ and the
amplitude is $E_{0}=1~{}{\rm V/m}$, while $h=\lambda/6$ and the incident angle
$\theta_{\rm i}=0\degree$.

The model of a conventional reflectarray corresponds to the case when
$d=\lambda/2$ (an array of half-wavelength spaced antenna elements loaded by
controllable loads). In this case, the solution for loads that ensure
excitation of the desired currents (8) and (11) is unique. Importantly, we
find that the required load impedances are complex-valued. That is, although
the array is overall lossless, individual antenna elements should either
receive or radiate power. This property is well known for anomalously
reflected periodical metasurfaces realized as impedance boundaries whose
impedance is defined by the desired plane-wave field distribution [10, 11]. As
to reflectarrays, we have found a discussion of this active-lossy response
problem in a classic book [35]. Obviously, the realizations of active and
lossy antenna elements in reflectarrays are not practically reasonable. For
periodical metasurfaces, the solutions have been found in ensuring the
excitation of optimized distributions of evanescent fields. Basically, these
evanescent fields form surface waves that effectively carry power from the
elements that need to receive power to the elements that need to radiate
power, in such a way that the perfect functionality with the purely reactive
local response of each array element is achieved. For reflectarrays, this
problem has not been even properly recognized, and no solutions are known.

In order to find the solutions for perfect anomalous reflection, first we need
to understand how significant is the effect of the required active or lossy
response of array elements. We consider the case where the angle of incidence
is zero and study the anomalous reflection efficiency as a function of the
desired reflection angle. For the calculated complex-valued loads the
efficiency obviously equals unity for all reflection angles, by the
construction of the solution. To see the effects of the active-lossy response,
we set the real parts of all calculated load impedances to zero and study how
the anomalous reflection efficiency changes. The corresponding induced
currents are found using (5), and the efficiency is obtained from (13).
Efficiency $\zeta$ as a function of the reflection angle $\theta_{\rm r}$ is
shown in Fig. 2 (see blue circles).

Figure 2: Efficiency as a function of reflection angle $\theta_{\rm r}$. The
red line is the theoretical limit for phase-gradient reflectors [12]. The
scatter circles are calculated at a step of $1\degree$ from $0\degree$ to
$89\degree$, when the real parts of the load impedances for $\lambda/2$-spaced
array are set to zero. The cross mark ’$\times$’ is the efficiency of the
conventional LPA method. The positive mark ’+’ represents the optimization
results that have only one strip in $\lambda/2$-spaced cells ($M=1$). The red
stars are the optimization results which have $3$ strips in $\lambda/2$-spaced
cell in Section IV when desired reflection angles are $65\degree$,
$70\degree$, and $75\degree$, respectively.

As a reference, we plot also the curve showing the fundamental limit for the
efficiency of phase-gradient reflectors defined by [12]

$\zeta=\frac{4\cos\theta_{\rm i}\cos\theta_{\rm r}}{(\cos\theta_{\rm
i}+\cos\theta_{\rm r})^{2}}$ (14)

(see red solid line in Fig. 2). For comparison, we calculate the required load
reactances in the conventional locally periodical approximation (LPA) and the
efficiency of the corresponding phase-gradient reflector. To do that, we
numerically calculate the reflection phase of a periodical array formed by one
loaded wire in a $\lambda/2$ cell at normal illumination. Next, we form a
phase-gradient reflector loading each wire so that a linearly varying profile
is realized in this locally periodical approximation. The efficiencies for
several reflection angles are shown by red crosses. It is known that for
phase-gradient metasurfaces with subwavelength cells, the efficiency of LPA
designs fall below the fundamental limit (14) (a collection of experimental
data from the literature can be found in [11]). We see that this conclusion
holds true also for reflectarrays with $\lambda/2$ unit cells.

Considering the results for analytically calculated load impedances with
dropped real parts, we observe that for angles $\theta_{\rm r}<30\degree$, the
efficiency of anomalous reflection is nearly ideal, and for reflection angles
less than $45\degree$, the efficiency remains greater than $95\%$. Thus, this
design significantly outperforms phase-gradient reflectors. When the desired
reflection angle is small (smaller than about $75\degree$), the efficiency
decreases monotonically. For larger reflection angles, efficiency decreases
dramatically. Due to the finite size of the strip array, the efficiency
fluctuates up and down when the reflection angle is larger than about
$75\degree$, at some point going beyond the fundamental limit for phase-
gradient boundaries.

To better demonstrate the effect of the real part of the loads, let us compare
two examples in detail: $\theta_{\rm r}=30\degree$ and $\theta_{\rm
r}=70\degree$.

###  III-A $\theta_{\rm r}=30\degree$

In this case, the induced current amplitudes determined by (8) and (11) equal
$I_{\alpha}=j4.5975\times 10^{-5}\text{A}$,
$\left|I_{\beta}\right|=4.7045\times 10^{-5}\text{A}$. In this case, the phase
of the induced current is rather special, it has a periodic profile since the
chosen period $d=\lambda/2$ matches with the period of periodical reflectors
operating at the angle $\theta_{\rm r}=30\degree$. The nearly periodic
behavior of load impedances and induced currents in the strips can be seen in
Fig. 3. It is worth noting that the currents in the strips remain almost
unchanged after removing the real part of the loads for such a moderate tilt
of the reflected wave (see Fig. 3(b)).

Figure 3: (a) Real and imaginary parts of the load impedance as the function
of the strip number when $\theta_{\rm r}=30\degree$. The solid red line
corresponds to the real part of the exact complex values of the impedances.
The subfigure gives detailed information about the load impedance densities in
the black dotted box. The solid blue line corresponds to purely imaginary load
impedances. (b) Amplitude and phase of the current distribution with and
without the real part of the load impedance densities. The blue line
represents the case of complex-valued load impedances, while the red line
represents the solution with dropped real parts. We can see that the current
has nearly the same distribution both in amplitude and phase. Figure 4:
Scattered field distribution $\Re\\{E_{x}^{\rm sca}\\}~{}[{\rm V/m}]$, where
(a) and (b) are ideal solutions with active-lossy loads, (c) and (d) show the
results when the real part of the loads are dropped, (e) and (f) are designed
with the LPA method. (a), (c), and (e) are for the case when $\theta_{\rm
r}=30\degree$, (b), (d), and (f) are for the case when $\theta_{\rm
r}=70\degree$, respectively.

Although the induced currents and the corresponding load impedances are
calculated analytically, it is not convenient to use the analytical
expressions for the scattered fields, because, for finite arrays over an
infinite ground plane, the scattered field contains the plane wave reflected
from the ground plane outside of the array area. For this reason, for the
optimized arrays we compute the fields also numerically, using COMSOL
software. In the numerical setup, we simulate the case of a finite-sized
ground plane, which allows us to clearly show the distribution of the
anomalously reflected fields in space. Details of the simulation setup are
given in Appendix B.

The real part of the scattered electric field $x$-component distribution is
depicted in Figs. 4 (a) and (c). It is clearly seen that the scattered fields
have nearly the same distribution for exact active-lossy load impedances after
dropping the real part of load impedances. This property allows the simple and
analytical design of anomalous reflectors based on the known impedance matrix
of the array. For the conventional LPA design method, it is necessary to
calculate the phase of the reflection coefficient for varying loads. In
addition, the scattered electric field contains parasitic propagating modes.
Using the proposed direct solution for required load impedances, one only
needs to use (5) and perform simple matrix operations to find the loads and
discard the real parts of the obtained values. Compared with the conventional
LPA method, this design method is simpler and more efficient (see Fig. 2).

###  III-B $\theta_{\rm r}=70\degree$

Figure 5: (a) Real and imaginary parts of the load impedance as the function
of the strip number when $\theta_{\rm r}=70\degree$. The solid red line
corresponds to the real part of the exact complex values of the impedances,
while the solid blue line corresponds to the purely imaginary load impedance
densities. (b) Amplitude and phase of the current distribution with and
without the real part of load impedance densities. The blue line represents
the case of complex-valued load impedances, while the red line represents the
solution with dropped real parts.

However, it is not always the case. For significantly larger reflection
angles, after directly discarding the real part of the loads, the induced
currents will be very different, resulting in a drop in efficiency, as part of
the energy will be scattered predominantly to the specular direction.

For example, when the desired reflection angle is $\theta_{\rm r}=70\degree$,
the induced currents determined through (8) and (11) are equal to
$I_{\alpha}=j4.5975\times 10^{-5}\text{A}$ and
$\left|I_{\beta}\right|=6.6424\times 10^{-5}\text{A}$. The load impedance
density distribution as a function of the strip number is depicted in Fig.
5(a) while the current distribution is depicted in Fig. 5(b). The induced
current distribution is quite different before and after neglecting the real
part of the load impedance densities. This obvious change in induced current
has a direct impact on the scattered field (see Figs. 4(b) and (d)).

The real part of the scattered electric field $x$-component distribution is
depicted in Figs. 4(b) and (d). We can see that when the real part of the load
impedance density is ignored, the reflected electric field is significantly
distorted and the efficiency drops. Constrained by power conservation, when
the real part of the load impedance density is removed, part of the energy is
scattered, mostly in the specular direction.

Although this design significantly outperforms conventional phase-gradient
reflectors, in this case the efficiency is significantly smaller than unity,
and more advanced approaches are needed, which we will introduce next.

##  IV Optimized subwavelength arrays to realize perfect anomalous reflectors

For large deviations from the usual reflection law (in the considered example
of reflection from normal incidence to a tilted angle larger than about
$50\degree$), the distribution of the induced current will be quite different
before and after neglecting the real part of the required load impedances. The
distorted current distribution deteriorates the radiation pattern in the far
zone due to the excitation of parasitic propagating modes. In this section, we
show that it is possible to dramatically improve the performance of the
reflector using optimization of current distribution over unit cells.

As already discussed, for periodical reflectors with the period determined by
the desired angles of incidence and reflection, the fundamental limit (14) can
be overcome by engineering distribution of near-zone (reactive) fields excited
by the incident wave [20, 36]. These known methods, unfortunately, are not
applicable for aperiodical, scanning reflectarrays that we consider here.
However, it is possible to engineer near-field distributions in Floquet-
periodical structures, where the current distribution along the array plane is
of the form

$I(y)=F(y)e^{-jk_{0}\sin\theta\,y}.$ (15)

Here, $F(y)$ is a periodical function with the period $\lambda/2$ (or
smaller). As is known from the theory of phased-array antennas (e.g., [33]),
such current distributions over infinite arrays launch a single plane wave
along the direction of the angle $\theta$. By optimizing function $F(y)$, it
is possible to control the distribution of evanescent fields in the vicinity
of the array, in analogy with the use of subwavelength meta-atoms in
periodical metasurfaces, but for scanning arrays with a fixed geometry of the
array.

In the frame of our simple model of arrays of thin impedance strips, to gain
control over the current distribution in each $\lambda/2$-spaced cell of the
array and be able to engineer the evanescent waves, we introduce additional
strips in each $\lambda/2$-spaced cell of the structure. It is assumed that
$M$ strips are equally spaced in each $\lambda/2$-spaced cell, so that the
distance between strips is $\lambda/(2M)$. This approach (with $M>1$) provides
the necessary freedom for control of current distribution over unit cells and
makes it possible to find loads with purely reactive impedances that can
realize perfect anomalous reflection in any direction. The averaged currents
in each $\lambda/2$ supercell are set to satisfy (8) and (11) (with
$d=\lambda/2$), so that the specular reflection is properly suppressed and the
power conservation in anomalous reflection is satisfied. However, having more
than one loaded strip in each unit cell, we can now optimize the subwavelength
distribution of the amplitude and phase of the induced current.

The induced current at the strip numbered $m$ can be expressed in the Floquet-
periodic form as

$I_{m}=F_{\alpha}(y_{m})I_{\alpha}e^{-jk_{0}\sin\theta_{\rm
i}D_{m}}+F_{\beta}(y_{m})I_{\beta}e^{j\phi}e^{-\text{j}k_{0}\sin\theta_{\rm
r}D_{m}},$ (16)

where $F_{\alpha}(y_{m})$ and $F_{\beta}(y_{m})$ are complex coefficients that
define how the total required currents $I_{\alpha}$ and $I_{\beta}$ in each
$\lambda/2$-supercell of the array are distributed between $M$ strips in the
supercell. We assume that the anomalously reflected wave has an arbitrary
phase shift with respect to the incident field, which we denote as $\phi$.
Functions $F_{\alpha}(y_{m})$ and $F_{\beta}(y_{m})$ are periodical with the
period $\lambda/2$. $D_{m}={\rm floor}(\frac{m}{M})\frac{\lambda}{2}$ defines
the linear phase gradient along the array of supercells. The function ${\rm
floor}(\cdot)$ returns the smallest integer not greater than its argument.

The constraints of the optimization are set on the averaged values over
supercells, to ensure power conservation. Due to the periodicity of
$F(\cdot)$, it is sufficient to give the constraints in one of the
$\lambda/2$-sized cells. Here, we give the constraints for the first
supercell. The constraints are nonlinear equality constraints and they follow
∑_m=0^M-1ℜ{F_α(y_m)} =& 1
∑_m=0^M-1ℑ{F_α(y_m)} = 0
∑_m=0^M-1ℜ{F_β(y_m)} = 1
∑_m=0^M-1ℑ{F_β(y_m)} = 0. These constraints ensure that the averaged currents
in each $\lambda/2$-sized cell have the same magnitude, while the phase of the
averaged currents in different $\lambda/2$-sized cells has a linearly varying
phase profile, as required for Floquet-periodical current distributions (15).

Most importantly, the periodicity of the distribution function $F(y)$ allows
us to avoid the need for global optimization of aperiodic arrays. It is enough
to optimize the current distribution over only one $\lambda/2$-sized unit
cell. Due to the periodicity of $F(y)$, the optimized distribution over one
cell determines the required currents induced in all strips of the whole array
(15). Finally, the corresponding load impedances for all strips are found
analytically using (5) based on the known impedance matrix of the array.

Thus, instead of more conventional direct optimizations of loads, here we
optimize the induced currents in each strip of one unit cell. Optimization is
performed by using function _fmincon_ in MATLAB. Our goal is to reroute all
the incident power in the desired direction. The objective function is to
maximize the amplitude of the far-zone electric field in the desired
direction. Once the goal is realized, the incident power is rerouted to the
desired direction in the most effective way, _i.e._ , perfect anomalous
reflection is realized.

The optimization process is organized as follows: First, an initial set of
currents subject to the nonlinear constraints is chosen. Then (5) is used to
calculate the corresponding load impedances. Since we optimize the induced
currents instead of the loads, the obtained loads are in general complex
numbers. With the goal to find an optimal solution with purely reactive loads,
we set the objective function as

$O=-{\rm max}\left\\{\left|E^{\rm far}(\theta_{\rm
r},\vec{I}(\Im\\{\vec{Z}_{\rm L}\\})\right|^{2}\right\\}.$ (17)

That is, we skip the real parts of the obtained load impedances and calculate
the corresponding induced currents. These current values are used to calculate
the scattered electric field in the far zone and the objective function. The
above steps are repeated until the termination conditions are met.

Let us discuss how many optimization parameters are needed. Because Matlab
cannot optimize complex numbers, we count the corresponding real-valued
parameters, the real and imaginary parts of the current amplitudes. For
periodic structures, the optimization is done in one period, requiring $2$
degrees of freedom for each propagating diffraction order [17]. For the
introduced aperiodic realization, the optimization is done in one
$\lambda/2$-sized cell. Each set of the induced currents requires $2M$
parameters since the induced current amplitudes are complex numbers. And we
have two sets of induced currents, one for eliminating the specularly
reflected wave, and another for generating the desired beam. Hence, the total
number of variables to be optimized is $4M+1$. An additional variable comes
from the phase of current $I_{\beta}$, which can be used in cases where the
application does not require control of the phase of the reflected beam.
Taking into account the nonlinear equality constraints, the actual number of
independent variables or degrees of freedom is $4M-3$. Compared with global
optimization, this optimization method has dramatically fewer variables to be
optimized, thus making the optimization speed and efficiency much higher,
especially when the total number of strips is quite large, for example, $N\geq
10$.

###  IV-A Optimized Subwavelength Strip Arrays

To verify the efficiency of optimization of current distribution inside
supercells, we give three examples of optimized anomalous reflectors. We
consider the case of $M=3$, that is, in every $\lambda/2$-spaced cell, we
position three impedance-loaded strips, so that the distance between the
strips is $d=\lambda/6$. There are $36$ supercells, each of which of
$\lambda/2$ width, correspondingly, the total width of the array is equal to
about $18\lambda$. The total number of strips equals $108$. In these examples,
we use the reflection phase as one of the optimization parameters, thus, we
optimize $4M-3=9$ real-valued variables. The array is illuminated by a
normally incident plane wave, and the desired reflection angles are
$65\degree,~{}70\degree$, and $75\degree$.

The optimization is done using the analytical expression for the scattered
fields (12), but for the graphical presentation we use numerical simulations,
utilizing the same setup as in Section III (Appendix B).

Figure 6: (a) - (c) Real part of the scattered electric field
$\Re\\{E_{x}^{\rm sca}\\}~{}[{\rm V/m}]$ distributions after optimization for
reflection angles $65\degree,70\degree$, and $75^{\degree}$, respectively. (d)
- (f) Distributions of the real part of the scattered electric field
$\Re\\{E_{x}^{\rm sca}\\}$ along the black solid arrow shown in subfigures
(a), (b), and (c). The red solid lines represent the theoretical optimal
amplitudes of the electric field for infinite arrays.

The optimized distributions of the scattered electric field $\Re{\\{E_{x}^{\rm
sca}\\}}$ for the three considered example reflection directions are shown in
Figs. 6 (a), (b), and (c), respectively. In Figs. 6 (d), (e), and (f), we also
plot the spatial field variations along the desired reflection direction,
shown by the arrows. We compare the field amplitudes with the theoretical
values of the required amplitudes of anomalously reflected plane waves from
perfect lossless infinite anomalous reflectors. These values are given by the
formula $E_{\rm r}=E_{\rm i}\sqrt{\cos\theta_{\rm i}/\cos\theta_{\rm r}}$ [8],
which returns $1.54,1.71$, and $1.97$ for the three considered reflection
angles. These values are shown by red lines on the corresponding plots. The
achieved far-field efficiencies, defined by (13), equal $99.4\%$, $109.3\%$,
$107.1\%$ for the considered example angles of reflection
$65\degree,~{}70\degree$ and $75\degree$. In these optimizations, the
reflection phase $\phi$ is considered as one of the optimization parameters.
For the given examples, the phase shifts of the optimized solutions are equal
to $\phi=~{}$0.0°$,~{}$255.5°$,~{}$43.8° for the angles
$65\degree,~{}70\degree$, and $75\degree$.

The results show excellent performance of the arrays of supercells with
optimized distributions of induced currents between three elements. The
incident power is rerouted to the desired direction, _i.e._ , perfect angle-
tunable anomalous reflection is realized using finite-sized arrays with a
fixed geometrical period. We also conclude that the optimization leads to
weakly super-directive arrays since the efficiency is above unity. In the last
example of reflection to $75\degree$ we see a weak sidelobe at a smaller
angle, yet, the far-field efficiency in the desired direction is above unity.
Interestingly, it appears that it is possible to realize nearly unit
efficiency for all angles even not using the reflection phase as a free
variable (fixing the desired phase), but this feature needs more studies.

###  IV-B Sparsely Arranged Strip Arrays

Finally, we compare the performance of the optimized subwavelength supercells
with arrays having only one controllable element per $\lambda/2$ cell. In
Section III, we already studied such arrays, analytically calculating the load
impedances that realize the desired distribution of induced currents. We see
that these impedances are complex-valued, and the performance after skipping
the real parts of the load impedances is better than for the conventional
phase-gradient reflector. However, this approach to finding the load
reactances does not ensure optimal performance, and it is interesting to study
if the performance of conventional arrays with a half-wavelength period can be
improved by using the proposed optimization approach.

Figure 7: (a) and (b) Scattered electric fields $\Re\\{E_{x}^{\rm
sca}\\}~{}[{\rm V/m}]$ for the case $\theta_{r}=55\degree$, when the real part
of the loads are neglected. (a) Analytical (non-optimized) solution, (b) after
optimization. (c) Induced current distribution. Red lines show data for
complex-valued load impedances and blue curves present the results for
reactive loads (real parts dropped).

For a sparsely distributed strip array $M=1$, _i.e._ , when there is only one
strip in $\lambda/2$-spaced cell. This means that only one degree of freedom
(the phase of $I_{\beta}$) can be adjusted to improve reflection efficiency.

The results of reflection efficiency after optimization are shown in Fig. 2
with the ’+‘ mark. We see that the reflection efficiency achieved with purely
reactive loads, in this case, depends on the phase of $I_{\beta}$, that is, on
the relative phase of the reflected field. If there is no need to control this
phase, the reflection efficiency can be further improved by properly selecting
the phase of $I_{\beta}$.

For the case $\theta_{\rm r}=55\degree$, the reflection efficiency is enhanced
from $77.3\%$ to $94.0\%$. It can be seen from Figs. 7(a) and (b), that the
scattered electric field is more “clean” after optimization which suggests a
higher efficiency. This result can be explained by considering the
corresponding induced current distribution which is depicted in Fig. 7(c).
After optimizing the phase of $I_{\beta}$, the difference between the case
with the complex-valued load impedances and the case with the neglected real
parts becomes smaller compared with the case without optimization both in the
amplitude and phase. Figure 7(c) shows that the current distribution after
optimization better resembles the ideal one compared to the current
distribution obtained by simply dropping the real parts of the load
impedances. As a result, the induced currents after optimization provide a
higher reflection efficiency.

We stress that for this case of a single strip in every $\lambda/2$ cell, we
cannot realize perfect anomalous reflection. The reason is that we cannot
engineer evanescent waves, which are necessary to match the impedances of the
incident and reflected waves. As discussed and illustrated above, this problem
can be overcome by using densely arranged strip arrays and other means to
control the current distribution at the subwavelength scale.

##  V Conclusions

In summary, we have systematically studied the means to realize finite-sized
reflectarrays capable of scanning the reflected beam to any angle without
parasitic scattering. The main motivation is the fact that the known solutions
for the realization of perfect anomalous reflectors are known only for
periodical arrays and metasurfaces that do not allow scanning of the
reflection angle. The study is performed using a simple, analytically solvable
two-dimensional model of an array of thin strips or wires. First, considering
an array of strips with a period $\lambda/2$ (that is, a model of a
conventional scanning reflectarray), we calculate analytically the load
impedances that need to be connected to each array element in order to realize
the desired distribution of induced currents (uniform amplitude and linear
phase gradient). The result shows that these loads must be active or lossy.
That is, for perfect performance the reflection must be nonlocal. At
effectively lossy cells, some part of the incident power is accepted by the
array but not re-radiated at the same point. Instead, this power must be
transported to another (effectively active) element and be re-radiated there.
This is the same phenomenon as is known for the realization of anomalous
reflection by using periodic metasurfaces. For scanning reflectarrays with
aperiodic distributions of load reactances, this issue was noticed in a
classic publication [35], but the means to overcome this limitation have not
been properly explored. It is important to note that this result also explains
the physical reason for parasitic specular reflection in off-set
reflectarrays.

One simple possibility that we have considered is to drop the real parts of
the required load impedances and use the calculated imaginary parts to set the
loads for each element. Naturally, the performance of the array is not anymore
perfect, but the results show that it is significantly better than the
conventional reflectarrays designed using the local periodic approximation.
For reflectors designed to reflect normally incident waves towards arbitrary
directions, the reflection efficiency is above $95\%$ for the reflection
angles smaller than $45\degree$. The required load reactances are found
algebraically, based on the pre-calculated impedance matrix of the array.

We show also that for a conventional array with $\lambda/2$ cells, the
performance of the reflector can be improved by optimizing the phase of the
reflected field, but, in that case, it is still impossible to achieve perfect
anomalous reflection at extreme angles, because of impedance mismatch between
the incident and reflected waves. For periodic metasurfaces, the problem of
the required nonlocal response is solved by optimizing surface waves (higher-
order Floquet harmonics). Here, we introduce and study the possibility to
achieve perfect anomalous reflection by optimizing current distribution over
$\lambda/2$-supercells that contain more than one strip. We have demonstrated
that in this way evanescent waves can be engineered for perfect anomalous
reflection at any angle. The simulation results have shown the feasibility and
effectiveness of this method. Importantly, no global optimization of load
reactances is needed, as it is enough to optimize the relative amplitudes and
phases of currents inside only one supercell of $\lambda/2$ size. Moreover,
the optimization is algebraic after the array impedance matrix is computed (in
this work, we used analytical expressions for mutual impedances). In the
presented example, perfect performance at all angles, even extreme, is
achieved by optimizing only 9 real-valued parameters for a linear aperiodic
array containing 108 elements.

Although the subwavelength-scale control is realized by using several
subwavelength antenna elements, it is possible to achieve the same goals by
using one single patch in each $\lambda/2$ unit cell by introducing several
feeding connectors to each patch, so that connecting different reactive loads
to different feeds, current distribution over the patch can be controlled.

Among other benefits, the proposed design method is not constrained to
illumination by a single plane wave and reflection towards one desired
direction. It can be also used to design reconfigurable beam splitters to
achieve different ratios of power distribution between different reflection
directions, and other devices for reflection control. Although we demonstrate
our method for TE polarized plane waves and scans in one plane only, the
physical principles are general, and this approach can be extended for TM
polarization and arbitrary scan directions in both planes. The main target
application is in reconfigurable intelligent surfaces for microwave,
millimeter, and THz ranges, but the same principles are applicable in optics.

##  Appendix A LPA method-based anomalous reflection

We calculate numerically the phase of the reflection coefficient from a
periodical array of identically loaded wires as a function of load reactance.
Using these results, we form a phase-gradient array loading the wires by
reactances needed to realize the desired local reflection coefficient phase.
Different loads are assigned to each strip according to different reflection
phase requirements at different positions. Next, we use the relationship
between the external electric fields at the wire positions, the load
impedances, and induced currents shown in (4) to calculate the currents
induced in all strips from

$\vec{I}=\overline{\overline{Z}}^{-1}\cdot\vec{U}.$ (18)

According to the known induced current distribution, we finally calculate the
scattered electric field and the anomalous reflection efficiency.

##  Appendix B Configuration of the aperiodic reflectarray numerical
simulation

The background field $\mathbf{E_{\rm b}}=\mathbf{E_{\rm
inc}}=e^{-jk_{0}\sin\theta_{\rm i}y-jk_{0}\cos\theta_{\rm i}z}\hat{x}$ is used
to excite the system. We calculate the corresponding loads using (5) and
insert these values into commercial software COMSOL Multiphysics to calculate
the scattered field. The configuration of COMSOL simulation setting is
displayed in Fig. 8.

Figure 8: Schematic of the 2D Comsol Multiphysics simulation configuration.
The blue area is the perfect matched layer (PML). The surface current density
on the wires is defined in terms of an impedance boundary condition, as
$\mathbf{J}=E_{x}/\left(Z_{{\rm L},n}/(2\pi r_{\rm eff})\right)\hat{x}$.

## References

  * [1] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. De Rosny, and S. Tretyakov, “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [2] A. Díaz-Rubio, S. Kosulnikov, and S. A. Tretyakov, “On the integration of reconfigurable intelligent surfaces in real-world environments: A convenient approach for estimation reflection and transmission,” _IEEE Antennas and Propagation Magazine_ , vol. 64, no. 4, pp. 85–95, 2022.
  * [3] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Communications Magazine_ , vol. 58, no. 1, pp. 106–112, 2019.
  * [4] D. Berry, R. Malech, and W. Kennedy, “The reflectarray antenna,” _IEEE Transactions on Antennas and Propagation_ , vol. 11, no. 6, pp. 645–651, 1963.
  * [5] S. V. Hum, M. Okoniewski, and R. J. Davies, “Modeling and design of electronically tunable reflectarrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 55, no. 8, pp. 2200–2210, 2007.
  * [6] J. Huang and J. A. Encinar, _Reflectarray Antennas_. Wiley, 2008.
  * [7] N. Yu, P. Genevet, M. A. Kats, F. Aieta, J.-P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: generalized laws of reflection and refraction,” _Science_ , vol. 334, no. 6054, pp. 333–337, 2011.
  * [8] V. S. Asadchy, M. Albooyeh, S. N. Tcvetkova, A. Díaz-Rubio, Y. Ra’di, and S. Tretyakov, “Perfect control of reflection and refraction using spatially dispersive metasurfaces,” _Physical Review B_ , vol. 94, no. 7, p. 075142, 2016.
  * [9] A. Epstein and G. V. Eleftheriades, “Synthesis of passive lossless metasurfaces using auxiliary fields for reflectionless beam splitting and perfect reflection,” _Physical Review Letters_ , vol. 117, no. 25, p. 256103, 2016.
  * [10] D.-H. Kwon, “Lossless scalar metasurfaces for anomalous reflection based on efficient surface field optimization,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 17, no. 7, pp. 1149–1152, 2018.
  * [11] A. Díaz-Rubio, V. S. Asadchy, A. Elsakka, and S. A. Tretyakov, “From the generalized reflection law to the realization of perfect anomalous reflectors,” _Science Advances_ , vol. 3, no. 8, p. e1602714, 2017.
  * [12] V. S. Asadchy, A. Wickberg, A. Díaz-Rubio, and M. Wegener, “Eliminating scattering loss in anomalously reflecting optical metasurfaces,” _ACS Photonics_ , vol. 4, no. 5, pp. 1264–1270, 2017.
  * [13] T. He, T. Liu, S. Xiao, Z. Wei, Z. Wang, L. Zhou, and X. Cheng, “Perfect anomalous reflectors at optical frequencies,” _Science Advances_ , vol. 8, no. 9, p. eabk3381, 2022.
  * [14] Y. Ra’di and A. Alù, “Reconfigurable metagratings,” _ACS Photonics_ , vol. 5, no. 5, pp. 1779–1785, 2018.
  * [15] A. Epstein, “Metagratings for perfect anomalous refraction,” in _2018 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting_. IEEE, 2018, pp. 923–924.
  * [16] V. Popov, F. Boust, and S. N. Burokur, “Beamforming with metagratings at microwave frequencies: Design procedure and experimental demonstration,” _IEEE Transactions on Antennas and Propagation_ , vol. 68, no. 3, pp. 1533–1541, 2019.
  * [17] V. Popov, F. Boust, and S. N. Burokur, “Constructing the near field and far field with reactive metagratings: Study on the degrees of freedom,” _Physical Review Applied_ , vol. 11, no. 2, p. 024074, 2019.
  * [18] V. Popov, M. Yakovleva, F. Boust, J.-L. Pelouard, F. Pardo, and S. N. Burokur, “Designing metagratings via local periodic approximation: From microwaves to infrared,” _Physical Review Applied_ , vol. 11, no. 4, p. 044054, 2019.
  * [19] A. Casolaro, A. Toscano, A. Alù, and F. Bilotti, “Dynamic beam steering with reconfigurable metagratings,” _IEEE Transactions on Antennas and Propagation_ , vol. 68, no. 3, pp. 1542–1552, 2019.
  * [20] X. Wang, A. Díaz-Rubio, and S. A. Tretyakov, “Independent control of multiple channels in metasurface devices,” _Physical Review Applied_ , vol. 14, no. 2, p. 024089, 2020.
  * [21] J. Budhu and Y. Rahmat-Samii, “Understanding the appearance of specular reflection in offset fed reflectarray antennas,” in _2011 IEEE International Symposium on Antennas and Propagation (APSURSI)_ , 2011, pp. 97–100.
  * [22] R. El hani and J.-J. Laurin, “Specular reflection analysis for off-specular reflectarray antennas,” _IEEE Transactions on Antennas and Propagation_ , vol. 61, no. 7, pp. 3575–3581, 2013.
  * [23] N. M. Estakhri and A. Alù, “Wave-front transformation with gradient metasurfaces,” _Physical Review X_ , vol. 6, no. 4, p. 041008, 2016.
  * [24] A. M. Wong and G. V. Eleftheriades, “Perfect anomalous reflection with a bipartite huygens’ metasurface,” _Physical Review X_ , vol. 8, no. 1, p. 011036, 2018.
  * [25] Y. Ra’di, D. L. Sounas, and A. Alù, “Metagratings: Beyond the limits of graded metasurfaces for wave front control,” _Physical Review Letters_ , vol. 119, no. 6, p. 067404, 2017.
  * [26] A. Epstein and O. Rabinovich, “Unveiling the properties of metagratings via a detailed analytical model for synthesis and analysis,” _Physical Review Applied_ , vol. 8, no. 5, p. 054037, 2017.
  * [27] O. Rabinovich and A. Epstein, “Analytical design of printed circuit board (pcb) metagratings for perfect anomalous reflection,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 8, pp. 4086–4095, 2018\.
  * [28] V. Popov, F. Boust, and S. N. Burokur, “Controlling diffraction patterns with metagratings,” _Physical Review Applied_ , vol. 10, no. 1, p. 011002, 2018\.
  * [29] O. Rabinovich, I. Kaplon, J. Reis, and A. Epstein, “Experimental demonstration and in-depth investigation of analytically designed anomalous reflection metagratings,” _Physical Review B_ , vol. 99, no. 12, p. 125101, 2019.
  * [30] G. Xu, G. V. Eleftheriades, and S. V. Hum, “Analysis and design of general printed circuit board metagratings with an equivalent circuit model approach,” _IEEE Transactions on Antennas and Propagation_ , vol. 69, no. 8, pp. 4657–4669, 2021.
  * [31] T. B. Hansen, “Non-uniform gratings that radiate preselected plane waves,” _IEEE Transactions on Antennas and Propagation_ , 2022.
  * [32] S. Tretyakov, _Analytical modeling in applied electromagnetics_. Artech House, 2003.
  * [33] F. Liu, D.-H. Kwon, and S. A. Tretyakov, “Reflectarrays and metasurface reflectors as diffraction gratings,” _arXiv preprint arXiv:2202.09029_ , 2022\.
  * [34] L. B. Felsen and N. Marcuvitz, _Radiation and scattering of waves_. John Wiley & Sons, 1994, vol. 31.
  * [35] O. Tereshin, V. Sedov, and A. Chaplin, _Synthesis of antennas on wave-delay structures (in Russian)_. Moscow: Svyaz, 1980.
  * [36] D.-H. Kwon, “Planar metasurface design for wide-angle refraction using interface field optimization,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 20, no. 4, pp. 428–432, 2021.