# Generalized Thermalization in an Integrable Lattice System

Amy C. Cassidy  Charles W. Clark  Joint Quantum Institute, National Institute
of Standards and Technology, Gaithersburg, MD 20899, USA  Marcos Rigol
Department of Physics, Georgetown University, Washington, DC 20057, USA

###### Abstract

After a quench, observables in an integrable system may not relax to the
standard thermal values, but can relax to the ones predicted by the
Generalized Gibbs Ensemble (GGE) [M. Rigol et al., PRL 98, 050405 (2007)]. The
GGE has been shown to accurately describe observables in various one-
dimensional integrable systems, but the origin of its success is not fully
understood. Here we introduce a microcanonical version of the GGE and provide
a justification of the GGE based on a generalized interpretation of the
eigenstate thermalization hypothesis, which was previously introduced to
explain thermalization of nonintegrable systems. We study relaxation after a
quench of one-dimensional hard-core bosons in an optical lattice. Exact
numerical calculations for up to 10 particles on 50 lattice sites ($\approx
10^{10}$ eigenstates) validate our approach.

###### pacs:

02.30.Ik,03.75.Kk,05.30.Jp,67.85.Hj

Once only of theoretical interest, integrable models of one-dimensional (1D)
quantum many-body systems can now be realized with ultracold atoms paredes04 ;
*kinoshita04. The possibility of controlling the effective dimensionality and
the degree of isolation have allowed access to the quasi-1D regime and to the
long coherence times necessary to realize integrable models. Additionally,
advances in the cooling and trapping of atoms have led to increased interest
in dynamics following quantum quenches, where a many-body system in
equilibrium is exposed to rapid changes in the confining potential or
interparticle interactions.

In general, in integrable quantum systems that are far from equilibrium,
observables cannot relax to the usual thermal state predictions because they
are constrained by the non-trivial set of conserved quantities that make the
system integrable sutherland04 . Relaxation to non-thermal values were
recently observed in a cold-atom system close to integrability kinoshita06 .
At integrability, it is natural to describe the observables after relaxation
by an updated statistical mechanical ensemble: the generalized Gibbs ensemble
(GGE) rigol07STATa , which is constructed by maximizing the entropy subject to
the integrability constraints jaynes57a ; *jaynes57b. In recent studies of
integrable systems rigol07STATa ; cazalilla06 ; *calabrese07a; *cramer08a;
*barthel08; *fioretto2010quantum; *rigol06STATb; kollar08 , the GGE has been
found to accurately describe various observables after relaxation, but a
microscopic understanding of its origin and applicability remains elusive. In
particular, an important question remains: how is it that expectation values
after relaxation can be described by an ensemble with exponentially fewer
parameters than the size of the Hilbert space? The full dynamics are
determined by as many parameters as the size of the latter. At a microscopic
level, thermalization for non-integrable systems can be understood in terms of
the eigenstate thermalization hypothesis (ETH) deutsch91 ; *srednicki94;
rigol08STATc , which, however, breaks down as one approaches integrability
rigol09STATa ; *rigol09STATb.

This paper is devoted to the study of how generalized thermalization, in the
sense of relaxation to the predictions of the GGE, takes place in integrable
systems. Answering this question is important not merely because of its
relevance to the foundations of statistical mechanics in integrable systems,
but also because it has become necessary to understand recent experiments with
ultracold gases in quasi-1D geometries. For integrable systems, we compare the
predictions of quantum mechanics with those of various statistical ensembles.
In particular, we introduce a microcanonical version of the GGE, which we use
to show that relaxation to the GGE can be understood in terms of a generalized
view of the ETH.

We study the dynamics following an instantaneous quench of 1D hard-core bosons
on a lattice, which is fully integrable. The Hamiltonian is given by

$\hat{H}=-J\sum_{i=1}^{L-1}\left(\hat{b}_{i}^{\dagger}\hat{b}_{i+1}+\textrm{H.c.}\right)+V(\tau)\sum_{i=1}^{L}(i-L/2)^{2}\hat{n}_{i}$
(1)

where $J$ is the hopping parameter; $V(\tau)$ gives the curvature of an
additional parabolic trapping potential for atoms on a lattice with lattice
constant $a$; $\hat{b}_{i}^{\dagger}\ (\hat{b}_{i})$ is the hard-core bosonic
creation (annihilation) operator; and
$\hat{n}_{i}=\hat{b}_{i}^{\dagger}\hat{b}_{i}$ is the number operator. In
addition to the standard commutation relations for bosons, hard-core bosons
satisfy the constraint $\hat{b}_{i}^{\dagger 2}=\hat{b}_{i}^{2}=0$, which
forbids multiple occupancy of the lattice sites. This Hamiltonian can be
mapped onto non-interacting fermions through the Jordan-Wigner transformation
jordan28 , and the many-body eigenstates can be constructed as Slater
determinants of the single-particle fermionic eigenstates rigol05HCBc .

We will focus on the behavior of the momentum distribution function,
$\langle\hat{n}_{k}\rangle=\sum_{l,m}e^{-ik(l-m)}\langle\psi|\hat{b}_{m}^{\dagger}\hat{b}_{l}|\psi\rangle/L$,
for system sizes ranging from $N=5$ bosons on $L=25$ lattice sites to $N=10$
bosons on $L=50$ lattice sites ($\approx 10^{10}$ eigenstates). Initially, we
prepare the system in the ground state $|\psi_{0}\rangle$ of a 1D lattice with
hard-wall boundary conditions and an additional harmonic potential, with
trapping strength $V=V_{0}$. At time $\tau=0$, the harmonic trap is turned
off, $V(\tau\geq 0)=0$, and the state $|\psi(\tau)\rangle$ evolves under the
influence of the final Hamiltonian. Hereafter, we refer to this state as it is
immediately after the quench as the “quenched state”. Its time evolution is
given by
$|\psi(\tau)\rangle=\sum_{\alpha}c_{\alpha}e^{-iE_{\alpha}\tau/\hbar}|\alpha\rangle,$
where $|\alpha\rangle$ are the energy eigenstates of the final Hamiltonian
with energies $E_{\alpha}$, and
$c_{\alpha}=\langle\alpha\vphantom{\psi_{0}}|\psi_{0}\vphantom{\alpha}\rangle$
are the overlaps between the eigenstates of the final Hamiltonian and the
quenched state. After relaxation, assuming the degeneracies in energy levels
are irrelevant, the expectation value of an observable is expected to be given
by the so called diagonal ensemble (DE) rigol08STATc ; kollar08 ; rigol09STATa
; *rigol09STATb

$\langle\hat{A}\rangle_{\text{DE}}=\lim_{\tau\rightarrow\infty}\frac{1}{\tau}\int_{0}^{\tau}d\tau^{\prime}\langle\psi(\tau^{\prime})|\hat{A}|\psi(\tau^{\prime})\rangle=\sum_{\alpha}|c_{\alpha}|^{2}\langle\alpha|\hat{A}|\alpha\rangle.$

We have checked numerically that, despite the integrability of our model,
$n_{k}$ relaxes to the DE prediction, with small fluctuations around this
result supp_mat .

Figure 1: (a),(b) Momentum distribution of the initial state (init), diagonal
(DE), generalized microcanonical (GME), generalized Gibbs (GGE), and the
microcanonical (ME) ensembles. (c),(d) Relative difference of the GME, GGE and
ME from the DE. (e),(f) Conserved quantities, $\langle\hat{I}_{n}\rangle$, in
the quenched state (identical to the DE), GME and ME. $\langle I_{n}\rangle$
are ordered in descending occupations in the quenched state. $L=50$, $N=10$,
$\delta_{\text{ME}}=0.05J$, $\delta_{\text{GME}}=0.8$. (a),(c),(e)
$\varepsilon=0.72J$, $V_{0}=0.029J$. (b),(d),(f) $\varepsilon=1.52J$,
$V_{0}=0.125J$.

Figure 1 shows the momentum distributions, $n_{k}$, before and after the
quench, for two different initial trap strengths, which correspond to
different energies per particle, $\varepsilon$, after the quench. These
results are compared with those of various ensembles of statistical mechanics.
The microcanonical ensemble (ME) is one in which all eigenstates in the
relevant energy window have identical weights. Within the microcanonical
ensemble, the expectation value of a generic observable $A$ is
$\langle\hat{A}\rangle_{\text{ME}}=N^{-1}_{\varepsilon,\delta_{\text{ME}}}\sum_{\alpha,|\varepsilon-\varepsilon_{\alpha}|<\delta_{\text{ME}}}\langle\alpha|\hat{A}|\alpha\rangle$,
where $\delta_{\text{ME}}$ is small, but still much greater than the mean
many-body level spacing. $N_{\varepsilon,\delta_{\text{ME}}}$ is the number of
eigenstates in the energy window
$|\varepsilon-\varepsilon_{\alpha}|<\delta_{\text{ME}}$. We have checked that
the results reported here are nearly independent of the specific value of
$\delta_{\text{ME}}$. The GGE is a grand-canonical ensemble that maximizes the
entropy subject to the constraints associated with non-trivial conserved
quantities of the quenched state. The density matrix takes the form
rigol07STATa

$\hat{\rho}_{\text{GGE}}=Z_{G}^{-1}e^{-\sum\lambda_{n}\hat{I}_{n}},\qquad
Z_{G}=\text{Tr}\left[e^{-\sum\lambda_{n}\hat{I}_{n}}\right],$ (2)

where {$\hat{I}_{n}$}, $n=1,\ldots,L$, are the conserved quantities. In our
systems, these correspond to the occupation of the single-particle eigenstates
of the underlying noninteracting fermions to which hard-core bosons are
mapped, and {$\lambda_{n}$} are Lagrange multipliers fixed by the initial
conditions,
$\lambda_{n}=\ln[(1-\langle\psi_{0}|\hat{I}_{n}|\psi_{0}\rangle)/\langle\psi_{0}|\hat{I}_{n}|\psi_{0}\rangle]$
rigol07STATa . Observables within this ensemble are then computed as
$\langle\hat{A}\rangle_{\text{GGE}}=\text{Tr}\left[\hat{A}\,\hat{\rho}_{\text{GGE}}\right]$
following Ref. rigol05HCBc .

As a step towards understanding the GGE as well as developing a more accurate
description of isolated integrable systems after relaxation, we introduce a
microcanonical version of the GGE, which we call the generalized
microcanonical ensemble (GME). Like the ME, where states within a small energy
window contribute with equal weight, within the GME we assign equal weight to
all eigenstates whose values of the conserved quantities are close to the
desired values. The expectation value of a generic observable within the
generalized microcanonical ensemble is given by
$\langle\hat{A}\rangle_{\text{GME}}=\mathcal{N}^{-1}_{\\{I_{n}\\},\delta_{\text{G
ME}}}\sum_{\alpha,\delta_{\alpha}<\delta_{\text{GME}}}\langle\alpha|\hat{A}|\alpha\rangle$,
where $\sum_{\alpha,\delta_{\alpha}<\delta_{\text{GME}}}$ is a sum over
eigenstates that are within the GME window and
$\mathcal{N}_{\\{I_{n}\\},\delta_{\text{GME}}}$ is the number of states within
that window and $\delta_{\alpha}$ is a measure of the distance of eigenstate
$\alpha$ from the target distribution.

In order to construct the GME, we include eigenstates of the Hamiltonian with
a similar distribution of conserved quantities which once averaged reproduce
the values of the conserved quantities in the quenched state. This approach is
characterized by three ingredients: (i) The ordered distribution (from largest
to smallest) of the conserved quantities in the DE, $\langle
I_{n}\rangle_{\textrm{DE}}\equiv\sum_{\alpha}|c_{\alpha}|^{2}I_{n,\alpha}$ [as
in Figs. 1(e) and 1(f)], (ii) a target distribution of the nonzero expectation
values of the conserved quantities $\\{I^{*}_{n^{*}_{i}}=1\\}$, where the
values of $n^{*}_{i}$ ($i=1,\ldots,N$) are chose to describe the distribution
$I_{n}$ in a coarse grained sense nstar , and (iii) for each individual many-
body eigenstate, the distance from the target state, $\delta_{\alpha}$, which
we define as
$\delta_{\alpha}=\left[\frac{1}{N}\sum_{i=1}^{N}I_{n^{*}_{i}}(n_{i,\alpha}-n^{*}_{i})^{2}\right]^{1/2}$.
Here $n_{i,\alpha}$ ($i=1,\ldots,N$) are the single-particle states occupied
in eigenstate $\alpha$, and $I_{n^{*}_{i}}$ are the interpolated values of
$\langle I_{n}\rangle_{\textrm{DE}}$, evaluated at $n^{*}_{i}$. The definition
of $\delta_{\alpha}$ is not unique and several variants that do not change our
conclusions were also considered supp_mat .

To better visualize the differences between the results of the various
ensembles in Figs. 1(a) and 1(b), we have plotted $\Delta\langle
n_{k}\rangle_{\text{stat}}=(\langle\hat{n}_{k}\rangle_{\text{DE}}-\langle\hat{n}_{k}\rangle_{\text{stat}})/\langle\hat{n}_{k}\rangle_{\text{DE}}$,
where “stat” stands for ME, GGE, or GME in Figs. 1(c) and 1(d). For weaker
initial confinements (smaller $\varepsilon$ \- Fig. 1(c)), the GME is
practically indistinguishable from the diagonal distribution. Both the GME and
the GGE accurately capture the tails of $n_{k}$, while the thermal ensemble
does not. For tighter initial traps (greater $\varepsilon$ \- Fig. 1(d)) all
four ensembles are very similar (note the scale), suggesting that $n_{k}$ in
the final steady state is indistinguishable from that of the thermal state.

The close agreement between DE and ME results in Fig. 1(b) raises the
question: how can an integrable system thermalize, given the constraints
imposed by the complete set of conserved quantities? We conjecture that if the
values of the conserved quantities in the quenched state are similar to those
of the ME, then the latter will accurately describe observables after
relaxation. This may occur for a variety of quenches.

In Figs. 1(e) and 1(f), we plot the values of the conserved quantities in the
quenched state and compare them with the expectation values of those
quantities in different statistical ensembles. (By definition, the
distribution of conserved quantities in the DE and GGE are identical to that
of the quenched state.) Figure 1(e) shows that the microcanonical values of
the conserved quantities are clearly different from the values in the quenched
state, while in Fig. 1(f) they are very similar. This supports our conjecture
above, and demonstrates that thermalization can occur in integrable systems
for special initial conditions. Additionally, the GME reproduces the correct
distribution of the conserved quantities supporting the validity of our method
for generating it.

To quantify the above observations, and to understand what happens in the
thermodynamic limit, we have studied the difference between the predictions of
the DE and the statistical ensembles for different system sizes. We compute
the integrated relative differences, $(\Delta
n_{k})_{\text{stat}}=\sum_{k}|\langle\hat{n}_{k}\rangle_{\text{DE}}-\langle\hat{n}_{k}\rangle_{\text{stat}}|/\sum_{k}\langle\hat{n}_{k}\rangle_{\text{DE}}$,
where again “stat” stands for ME, GGE, or GME.

Figure 2:  (a) $(\Delta n_{k})_{\text{ME}}$ versus energy per particle of the
quenched state. $\delta_{\text{ME}}=0.05J$. (b) $(\Delta n_{k})_{\text{GME}}$
vs $\varepsilon$, $\delta_{\text{GME}}=0.8$. (c) Integrated difference between
the conserved quantities in the quenched state and the ME, $(\Delta
I_{n})_{\text{ME}}$. (d) $(\Delta n_{k})_{\text{GGE}}$ vs $\varepsilon$.
Inset: $(\Delta n_{k})_{\text{GGE}}$ vs $L^{-0.73}$ for $\varepsilon=1.07J$,
where a fit to $(\Delta n_{k})_{\text{GGE}}=zL^{-\gamma}$ gives
$\gamma=0.73\pm 0.02$.

In Fig. 2(a), we plot $(\Delta n_{k})_{\text{ME}}$ as a function of the final
energy per particle, $\varepsilon$, for different lattice sizes, $L$. To
perform finite-size scaling, $\varepsilon$ and the filling factor
($\nu=N/L=0.2$) are held constant as $L$ changes. Figure 2(a) shows that for
$\varepsilon\lesssim 1.3J$ the difference between the $n_{k}$ in the DE and
the ME increases with increasing $L$, indicating that the difference persists
in the thermodynamic limit. For $\varepsilon\gtrsim 1.3J$, the opposite
behavior is observed. From our previous discussion, one expects that $(\Delta
n_{k})_{\text{ME}}$ should closely follow the behavior of the integrated
differences between the conserved quantities in the quenched state and the ME,
$(\Delta I_{n})_{\text{ME}}={\sum_{n}|\langle I_{n}\rangle_{\text{DE}}-\langle
I_{n}\rangle_{\text{ME}}|}/{\sum_{n}\langle I_{n}\rangle_{\text{DE}}}$. This
is seen by comparing Figs. 2(a) and 2(c), which leads us to conclude that
$n_{k}$ need not relax to the standard thermal prediction, except when
$(\Delta I_{n})_{\text{ME}}$ becomes negligible. Qualitatively similar results
were obtained in the canonical ensemble supp_mat .

On the other hand, in Fig. 2(b) one can see that the differences between
$n_{k}$ in the diagonal and generalized microcanonical ensembles are very
small and decrease with increasing system size, so that the former
successfully describes this observable after relaxation. In the case of the
GGE [Fig. 2(d)], $(\Delta n_{k})_{\text{GGE}}$ is in general larger than
$(\Delta n_{k})_{\text{GME}}$, which is to be expected since the GGE is a
grand-canonical ensemble. As the system size increases $(\Delta
n_{k})_{\text{GGE}}\rightarrow 0$ as $L^{-\gamma}$, where $\gamma\approx 0.73$
[inset of Fig. 2(d)] and slightly depends on the energy supp_mat .

The question that remains to be answered is why the generalized Gibbs and the
generalized microcanonical ensemble are able to describe the $n_{k}$ after
relaxation, i.e., why
$\langle\hat{n}_{k}\rangle_{\text{GGE}}=\langle\hat{n}_{k}\rangle_{\text{GME}}=\langle\hat{n}_{k}\rangle_{\text{DE}}\equiv\sum_{\alpha}|c_{\alpha}|^{2}\langle\alpha|\hat{n}_{k}|\alpha\rangle$.
Note that whereas $\langle\hat{n}_{k}\rangle_{\text{GGE}}$ and
$\langle\hat{n}_{k}\rangle_{\text{GME}}$ are entirely determined by the $L$
independent values of the conserved quantities in the quenched state,
$\langle\hat{n}_{k}\rangle_{\text{DE}}$ is determined by the exponentially
larger $\binom{L}{N}$ values of the coefficients $c_{\alpha}$.

To address this question, we perform a spectral decomposition of
$\langle\hat{n}_{k}\rangle_{\text{DE}}$ and
$\langle\hat{n}_{k}\rangle_{\text{GME}}$. Figure 3 displays a coarse grained
view of the weight which eigenstates with a given zero momentum occupancy
$\langle\hat{n}_{k=0}\rangle_{\alpha}=\langle\alpha|\hat{n}_{k=0}|\alpha\rangle$
contribute to the DE [Fig. 3(a)] and the GME [Fig. 3(b)]. The correlation
between the results in both figures is apparent. However, it is not clear why
the details contained in the overlaps $c_{\alpha}$ are completely washed out
so that the DE and the GME results coincide, while they are different from
those in the ME. In the inset of Fig. 3(a), we plot a histogram of the values
of $n_{k=0}$ for the DE, GME and the ME. Clearly the histograms for the DE and
GME have a similar mean but different widths, while the ME has a different
mean and width supp_mat .

Figure 3:  Density plot of the coarse-grained weights with which eigenstates
contribute to (a) the DE (sum over diagonal weights, $|c_{\alpha}|^{2})$ and
(b) the GME (fractional number of states = number of states/total number of
states) as a function of eigenstate energy and
$\langle\hat{n}_{k=0}\rangle_{\alpha}$. The sums are performed over window of
width $\delta n_{k}=0.0067$ , and $\delta\varepsilon=0.035J$. The horizontal
and vertical dotted lines are the expectation values of $\hat{n}_{k=0}$ and
$\varepsilon$ in each ensemble. $L=45$, $N=9$, $V_{0}=0.036J$,
$\varepsilon=0.72J$, $\delta_{\text{GME}}=0.85$. Inset in (a): Histogram of DE
weights (green), fractional number of GME states (blue) and fractional number
of ME states (black) summed over all energies. Bin width, $\delta
n_{k}=0.0067$. Vertical lines give the mean, $\langle\hat{n}_{k=0}\rangle$
within each ensemble. Inset in (b): Fluctuations of
$\langle\hat{n}_{ka=0}\rangle(\bullet)$ and
$\langle\hat{n}_{ka=2\pi/5}\rangle(\blacktriangle)$ within the DE (green) and
GME (blue) as a function of inverse system size. $\varepsilon=0.72J$.

Ultimately, one is interested in what happens in the thermodynamic limit. For
each $k$, we define the width of the distribution of
$\langle\hat{n}_{k}\rangle_{\alpha}$ for each ensemble as
$\sigma_{k}=\sqrt{\langle\hat{n}_{k}^{2}\rangle-\langle\hat{n}_{k}\rangle^{2}}$.
The inset of Fig. 3(b), shows $\sigma_{k}$ within the DE and the GME versus
$L^{-1}$. The scaling is depicted for two $k$ values and clearly shows that
the widths of both distributions vanish in the thermodynamic limit. This
demonstrates that the overwhelming majority of the states selected by the DE
as well as by the GME, which have similar values of the conserved quantities,
have identical expectation values of $n_{k}$. This is why details of the
distribution of $c_{\alpha}$ no longer matter as $L$ increases. We note that
with increasing $L$, the number of eigenstates contained in the generalized
microcanonical window increases exponentially, however, the ratio of the
number of states in the GME and the ME vanishes supp_mat .

The findings above provide a generalization of the ETH introduced previously
to understand thermalization in nonintegrable systems deutsch91 ;
*srednicki94; rigol08STATc . The ETH states that the expectation values of
few-body observables in generic systems do not fluctuate between eigenstates
that are close in energy. Thus all eigenstates within a microcanonical window
have essentially the same expectation values of the observables, and one can
say that thermalization occurs at the level of eigenstates. As seen in Fig. 3,
$\langle\hat{n}_{k}\rangle_{\alpha}$ exhibits large eigenstate-to-eigenstate
fluctuations in our integrable system, showing that ETH is invalid. However,
by selecting eigenstates with similar conserved quantities, ETH is restored,
although in a weaker sense: the overwhelming majority of eigenstates with
similar conserved quantities have similar values of $n_{k}$. These results
pave the way to a unified understanding of thermalization in generic
(nonintegrable systems) and its generalization in integrable systems. This
opens many new questions, such as whether the concepts of typicality tasaki98
; *goldstein06; *popescu06; *reimann08 and thermodynamics
polkovnikov08microscopic ; *polkovnikov08 can be generalized to isolated
integrable systems.

###### Acknowledgements.

This work was supported by NSF under Physics Frontier Grant No. PHY-0822671.
A.C.C. acknowledges support from NRC/NIST and M.R. acknowledges support from
the Office of Naval Research. We thank V. Dunjko, L. Mathey, and M. Olshanii
for helpful discussions.

## References

  * (1) B. Paredes _et al._ , Nature 429, 277 (2004)
  * (2) T. Kinoshita, T. Wenger, and D. Weiss, Science 305, 1125 (2004)
  * (3) B. Sutherland,  _Beautiful Models_ (World Scientific, Singapore, 2004)
  * (4) T. Kinoshita, T. Wenger, and D. S. Weiss, Nature 440, 900 (2006)
  * (5) M. Rigol _et al._ , Phys. Rev. Lett. 98, 050405 (2007)
  * (6) E. T. Jaynes, Phys. Rev. 106, 620 (1957)
  * (7) E. T. Jaynes, Phys. Rev. 108, 171 (1957)
  * (8) M. A. Cazalilla, Phys. Rev. Lett. 97, 156403 (2006)
  * (9) P. Calabrese and J. Cardy, J. Stat. Mech., P06008(2007)
  * (10) M. Cramer _et al._ , Phys. Rev. Lett. 100, 030602 (2008)
  * (11) T. Barthel and U. Schollwöck, Phys. Rev. Lett. 100, 100601 (2008)
  * (12) D. Fioretto and G. Mussardo, New J. Phys. 12, 055015 (2010)
  * (13) M. Rigol, A. Muramatsu, and M. Olshanii, Phys. Rev. A 74, 053616 (2006)
  * (14) M. Kollar and M. Eckstein, Phys. Rev. A 78, 013626 (2008)
  * (15) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991)
  * (16) M. Srednicki, Phys. Rev. E 50, 888 (1994)
  * (17) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008)
  * (18) M. Rigol, Phys. Rev. Lett. 103, 100403 (2009)
  * (19) M. Rigol, Phys. Rev. A 80, 053607 (2009)
  * (20) P. Jordan and E. Wigner, Z. Phys. 47, 631 (1928)
  * (21) M. Rigol, Phys. Rev. A 72, 063607 (2005)
  * (22) See supplementary materials online.
  * (23) In each many-body eigenstate of the Hamiltonian, $|\alpha\rangle$, $I_{n,\alpha}=\langle\alpha|\hat{I}_{n}|\alpha\rangle$ can only be zero or one, because they are constructed by occupying single-particle fermionic states.
  * (24) We compute the $n^{*}_{i}$, which are not restricted to integers, as follows. The first one ($n^{*}_{1}$) is computed so that $\int_{0.5}^{n^{*}_{1}}I(x)dx=0.5$, where $I(x)=I_{n},x\in(n-0.5,n+0.5]$. For all other values of $n^{*}_{i}$, $\int_{n^{*}_{i-1}}^{n^{*}_{i}}I(x)dx=1$.
  * (25) H. Tasaki, Phys. Rev. Lett. 80, 1373 (1998)
  * (26) S. Goldstein _et al._ , Phys. Rev. Lett. 96, 050403 (2006)
  * (27) S. Popescu, A. J. Short, and A. Winter, Nature Phys. 2, 754 (2006)
  * (28) P. Reimann, Phys. Rev. Lett. 101, 190403 (2008)
  * (29) A. Polkovnikov, Phys. Rev. Lett. 101, 220402 (2008)
  * (30) A. Polkovnikov, arXiv:0806.2862 (2008)

Supplementary material for EPAPS
Generalized Thermalization in an Integrable Lattice System

Amy C. Cassidy,1 Charles W. Clark,1 Marcos Rigol2

1Joint Quantum Institute, National Institute of Standards and Technology,
Gaithersburg, MD 20899, USA

2Department of Physics, Georgetown University, Washington, DC 20057, USA

Time Dynamics. In order to check that the DE accurately describes observables
after relaxation, we compare the time dynamics of the central momentum peak,
$\langle\hat{n}_{k=0}(\tau)\rangle=\langle\psi(\tau)|\hat{n}_{k=0}|\psi(\tau)\rangle$,
with the expectation value of $\hat{n}_{k=0}$ in the diagonal ensemble. Figure
1(a) shows that $\langle\hat{n}_{k=0}(\tau)\rangle$ relaxes to the diagonal
prediction, with small fluctuations around this result, indicating that the
diagonal ensemble correctly predicts the values of observables after
relaxation. Additionally, the expectation value of $\hat{n}_{k=0}$ in the GME
agrees with the DE, while the ME does not.

Figure 1: (a) Time evolution of $n_{k=0}$. Horizontal lines represent
$\langle\hat{n}_{k}\rangle_{\text{DE}}$,
$\langle\hat{n}_{k}\rangle_{\text{GME}}$,
$\langle\hat{n}_{k}\rangle_{\text{ME}}$. $L=50$, $N=10$, $\varepsilon=0.72J$,
$\delta_{\text{GME}}=0.8$, $\delta_{\text{ME}}=0.05J$. Main Panel: $\tau\leq
100$. Inset: $100\leq\tau\leq 1000$. (b) Integrated difference of the diagonal
and time-averaged momentum distribution,$(\Delta n_{k})_{\tau}$, versus
$L^{-0.98}$, where a fit to $\sigma_{\tau}=zL^{-\gamma}$ gives $\gamma=0.98\pm
0.25$. (c) Mean fluctuations of the momentum distribution, $\sigma_{\tau}$,
versus $L^{-0.4}$, where a fit to $(\Delta n_{k})_{\tau}=zL^{-\gamma}$ gives
$\gamma=0.40\pm 0.12$. (b),(c)
$\tau_{1}=100J,\tau_{2}=1000J,\varepsilon=0.72J$. $z$ is a generic
multiplicative constant.

In order to understand the relation between the DE prediction and the actual
time average over a finite time interval,
$\overline{n_{k}}=\frac{1}{\tau_{2}-\tau_{1}}\int_{\tau_{1}}^{\tau_{2}}d\tau\langle\hat{n}_{k}(\tau)\rangle$,
we study the integrated difference $(\Delta
n_{k})_{\tau}=\sum_{k}|\langle\hat{n}_{k}\rangle_{\text{DE}}-\overline{n_{k}}|/\sum_{k}\langle\hat{n}_{k}\rangle_{\text{DE}}$,
and the time fluctuations
$\sigma_{\tau}=\sum_{k}\sqrt{\overline{n_{k}^{2}}-\overline{n_{k}}^{2}}$ as a
function of system size. $(\Delta n_{k})_{\tau}$ is depicted in Fig. 1(b),
where we also show the results of a fit to $(\Delta
n_{k})_{\tau}=zL^{-\gamma}$ with $\gamma=0.98\pm 0.25$ ($z$ is a generic
multiplicative constant). The integrated difference, $(\Delta n_{k})_{\tau}$,
clearly decreases with increasing system size and is several orders of
magnitude smaller than similar comparisons with the statistical ensembles,
confirming that degeneracies are irrelevant and the diagonal distribution
accurately represents the long-time average. Figure 1(c) depicts the time
fluctuations of the momentum distribution, $\sigma_{\tau}$, where the exponent
was determined from fitting $\sigma_{\tau}=zL^{-\gamma}$. This plots makes
evident that the fluctuations also decrease with increasing system size
(likely $\propto 1/\sqrt{L}$) and are expected to vanish in the thermodynamic
limit.

Figure 2:  (a) Distribution of conserved quantities in the diagonal ensemble,
$\langle I_{n}\rangle_{\textrm{DE}}$, ordered by energy of the single-particle
fermionic eigenstates. (b) $\langle I_{n}\rangle_{\textrm{DE}}$ sorted in
descending order. The target distribution of $n_{i}^{\ast}$ are labeled as
well as the corresponding $I_{n_{i}^{\ast}}$ for $i=1,2$. (c) Target
distribution of conserved quantities, $I_{n}^{\ast}$, and distribution of
eigenstate that minimized the distance $\delta_{\alpha}$ used to generate the
GME, $I_{n,\alpha}$. (d) $\langle I_{n}\rangle$ in the DE and GME using method
described in text as well as the standard metric (GME-std).
$L=20,N=4,\varepsilon=0.72J$.

Constructing the generalized microcanonical ensemble. We expand upon our
method for generating the generalized microcanonical ensemble with a series of
plots in Fig. 2. Constructing the generalized microcanonical ensemble for a
system of 1D hard-core bosons presents some unique challenges. In particular,
the set of conserved quantities are the single-particle occupations of the
underlying fermions used to build the many body eigenstates. Thus in any given
eigenstate, each conserved quantity is either $0$ or $1$ corresponding to
whether or not that fermionic state is occupied, such as $I_{n,\alpha}$ in
Fig. 2(c). On the other hand, the distribution in the diagonal ensemble is a
continuous variable, $0\leq\langle I_{n}\rangle_{\textrm{DE}}\leq 1$. For each
individual eigenstate, we must include or exclude it from the GME based on how
close the discrete distribution in the eigenstate is to the continuous
distribution of the quenched state. One possible method is a simple extension
of the microcanonical metric

$\delta^{\prime}_{\alpha}=\left[\sum_{n=1}^{L}(I_{n}^{\ast}-I_{n,\alpha})^{2}/\sum_{n}I_{n,\alpha}^{2}\right]^{1/2}.$

We found this was not guaranteed to accurately describe the distribution of
conserved quantities, particularly for low energies for the system sizes
studied (see Fig. 2(d)). Instead of using the standard metric, we perform a
weighted least-squares fit to a target distribution which is described as
follows.

The unsorted conserved quantities of the quenched state $\langle
I_{n}\rangle_{\textrm{DE}}$, ordered by energy of the single-particle
fermionic eigenstates, are plotted in Fig. 2(a). As a first step, the
distribution of $\langle I_{n}\rangle_{\textrm{DE}}$ is sorted in descending
order as shown in Fig 2(b). Next the values of $n_{i}^{\ast}$, which are not
limited to integer values, are chosen so that
$\int_{0.5}^{n_{1}^{\ast}}I(x)=0.5$, where $I(x)=\langle
I_{n}\rangle_{\text{DE}}$ for $x$ in the interval $(n-0.5,n+0.5]$.
Subsequently the $n_{i}^{\ast}$ are determined so
$\int_{n_{i-1}^{\ast}}^{n_{i}^{\ast}}I(x)=1$ as depicted in Fig. 2(b). This
set of $\\{n_{i}^{\ast}\\}$ determine the target distribution, labeled
$I_{n}^{\ast}$ in Fig. 2(c). Each $n_{i}^{\ast}$ is assigned a weight,
$I_{n_{i}^{\ast}}$, which is the interpolated value of $\langle
I_{n}\rangle_{\textrm{DE}}$ at $n_{i}^{\ast}$. The distance between the
distribution of conserved quantities in eigenstate $\alpha$ and the
distribution of the quenched state is then defined as

$\delta_{\alpha}=\left[\frac{1}{N}\sum_{i=1}^{N}I_{n^{*}_{i}}(n_{i,\alpha}-n^{*}_{i})^{2}\right]^{1/2},$
(1)

where the $n_{i,\alpha}$ are the single particle eigenstates that are occupied
in the many-body eigenstate $\alpha$.

Our method of constructing the GME is not unique. We tested various different
approaches, including an unweighted sum in Eq. (1), etc, which gave similar
results. We then employed the method which best reconstructed the distribution
of conserved quantities over the full range of parameters studied. Note that
for clarity Fig. 2 displays data for the case of $N=4$ particles on $L=20$
sites. There are only 4845 total eigenstates and 174 eigenstates used to
construct the GME for the data shown. Even for such a small system, the
distribution of conserved quantities within the GME is quite good.

## Lagrange multipliers and additivity.

In Fig. 3, we plot the expectation values of the conserved quantities in the
diagonal ensemble along with corresponding Lagrange multipliers. The Lagrange
multipliers are given by $\lambda_{n}=\ln\left(\frac{1-\langle
I_{n}\rangle_{\textrm{DE}}}{\langle I_{n}\rangle_{\textrm{DE}}}\right)$. As
can be seen, the Lagrange multipliers vary smoothly with the value of $I_{n}$
and with the index $n$ after the conserved quantities have been ordered in
descending order. Additionally, the distribution of conserved quantities is
very similar for two different lattice sizes with the same final energy-per-
particle, when the index is normalized by the total number of conserved
quantities. This is an important property of the Lagrange multipliers, which
shows that even though the conserved quantities in these integrable systems
are not additive in a strict sense, they can still be understood to be
additive in a coarse grained sense, because the values of $\lambda_{n}$ are a
smooth function of $I_{n}$.

Figure 3: Expectation values of the conserved quantities in the diagonal
ensemble and corresponding Lagrange multipliers for $N=5$, $L=25$ and $N=10$,
$L=50$. $\varepsilon=0.72t$.

Canonical Ensemble. In addition to the results reported for the microcanonical
distribution, we study the momentum distribution in the canonical ensemble
(CE). Within the canonical ensemble, observables are calculated as
$\langle\hat{A}\rangle_{\text{CE}}=Z_{\text{CE}}^{-1}\text{Tr}\left[\hat{A}e^{-\beta\hat{H}}\right]$,
where $Z_{\text{CE}}=\text{Tr}\left[e^{-\beta\hat{H}}\right]$ and $\beta$ is
the inverse temperature. The temperature is calculated numerically so that
$\langle E\rangle_{\text{CE}}=\langle\psi_{0}|\hat{H}(\tau\geq
0)|\psi_{0}\rangle$. In Fig. 4(a) we plot the integrated difference between
the momentum distribution in the DE and CE, $(\Delta n_{k})_{\text{CE}}$, for
different lattice sizes as a function of the energy per particle of the
quenched state. The comparison of the canonical and diagonal momentum
distributions are similar to the comparison between the microcanonical and
diagonal distributions, although there is some discrepancy between the
microcanonical and canonical distributions due to finite size effects. In
particular, the upturn in $(\Delta n_{k})_{\text{ME}}$ at $\varepsilon=1.4J$
is not present in $(\Delta n_{k})_{\text{CE}}$.

Figure 4: (a) Integrated difference between the diagonal and canonical
momentum distributions, $\left(\Delta n_{k}\right)_{\text{CE}}$, versus energy
per particle, $\varepsilon$. (b) Integrated difference of the momentum in the
microcanonical and canonical ensembles, $\left(\Delta n_{k}\right)_{\text{CE-
ME}}$, versus $L^{-1.19}$ for $\varepsilon=1.07J$, where a fit to
$\left(\Delta n_{k}\right)_{\text{CE-ME}}=zL^{-\gamma}$ gives $\gamma=1.0\pm
0.07$. Inset: Scaling exponent, $\gamma$ of $\left(\Delta
n_{k}\right)_{\text{CE-ME}}$ as a function of the energy per particle,
$\varepsilon$.

In general, the microcanonical distribution provides us with a better
description of the system because it is isolated, although it breaks down for
very small systems because of poor statistics due to an insufficient number of
states in the relevant energy window. Given the expectation that
$\langle\hat{n}_{k}\rangle_{\text{CE}}$ and
$\langle\hat{n}_{k}\rangle_{\text{ME}}$ are equal in the thermodynamic limit,
we calculate the scaling exponent of the integrated difference of the momentum
distribution, $\left(\Delta n_{k}\right)_{\text{CE-ME}}=zL^{-\gamma}$, where
$\left(\Delta n_{k}\right)_{\text{CE-
ME}}=\sum_{k}|\langle\hat{n}_{k}\rangle_{\text{CE}}-\langle\hat{n}_{k}\rangle_{\text{ME}}|/\sum_{k}\langle\hat{n}_{k}\rangle_{\text{CE}}$.
In Fig. 4(b) we plot $\left(\Delta n_{k}\right)_{\text{CE-ME}}$ versus
$L^{-\gamma}$ for energy $\varepsilon=1.31J$, where $\gamma=1.0$. In the
inset, the exponent $\gamma$ is plotted as a function of the energy per
particle. The mean value is $\gamma=1.03\pm 0.04$, which is consistent with
$L^{-1}$ scaling for the difference between the canonical and grand canonical
momentum distributions for hard-core bosons in a box found in Ref. rigol05HCBc
.

On general grounds, the canonical and microcanonical ensembles are expected to
equivalent in the thermodynamic limit, which is confirmed by our numerical
results. We also see increasing agreement between the results of the
generalized Gibbs and generalized microcanonical ensembles as system sizes
increase. There has been significant work in recent years on the equivalence
of classical ensembles and less in quantum systems. In the classical case,
Ellis et al. have demonstrated that the microcanonical and canonical ensembles
are equivalent if and only if the thermodynamic functions are equivalent,
which is the case if the microcanonical entropy is concave ellis2000large .
Furthermore, when the two ensemble are not equivalent, it may be possible to
construct a generalized canonical ensemble, which contains an additional
exponential which is a continuous function of the Hamiltonian, that is
equivalent to the microcanonical ensemble costeniuc2005generalized .

Scaling of the GGE. Given the trend observed in our results, and similar
calculations in rigol07STATa  for much larger system sizes, we expect the
integrated difference between the momentum distribution in the DE and GGE,
$(\Delta n_{k})_{\text{GGE}}$, to vanish in the thermodynamic limit. We fitted
$(\Delta n_{k})_{\text{GGE}}=zL^{-\gamma}$ and report the results for
$\varepsilon=1.13J$ in the main text. In Fig. 5 we plot the exponent,
$\gamma$, as a function of the energy per particle of the quenched state for
all energies studied. The mean value of the exponent is $\gamma=0.731\pm
0.007$.

Figure 5: Exponent of power-law scaling of $(\Delta n_{k})_{\text{GGE}}$
versus inverse system size as a function of the energy per particle.

Number of States. The total number of states in our system, the number of
states in the ME, and the number of states in the GME all scale exponentially
with the system size as $z^{L}$. For the total number of states,
$L!/\left[N!(L-N)!\right]$, using Stirling’s approximation,
$z\approx(1-\nu)^{\nu-1}\nu^{-\nu}\approx 1.65$, where $\nu$ is the filling
factor. We confirm this scaling by numerically fitting to the total number of
states. For $\varepsilon=0.72J$, we find numerically that $z\approx 1.51$ in
the ME, ($\delta_{ME}=0.05J$) and $z\approx 1.37$ in the GME
($\delta_{GME}=0.8$). Thus the number of states in the ME and GME windows as a
the fraction of the total number of states vanishes in the thermodynamic
limit. Additionally the ratio of states in GME to ME vanishes. These values
are typical, although the precise exponents depend on the choice of the ME and
GME window sizes.

Figure 6: (a) Average density fluctuations, $\sigma_{\text{den}}$ versus
$L^{-0.5}$ for fixed energy per particle in the diagonal (DE), generalized
microcanonical (GME) and microcanonical (ME) ensembles. (b) Momentum
fluctuations, $\sigma_{k}$ for $ka=0$. $\varepsilon=0.72J$,
$\delta_{\text{ME}}=0.05J,\delta_{\text{GME}}=0.8.$

Fluctuations of Local Observables. For local observables, the fluctuation of
the eigenstate to eigenstate expectation values are expected to scale as
$L^{-1/2}$ biroli2009 . In Fig. 6(a) we plot the fluctuations of the site
occupations averaged over five sites,
$\sigma_{\text{den}}=1/5\sum_{i=1}^{5}\sqrt{\langle\hat{n}_{x=iN}^{2}\rangle-\langle\hat{n}_{x=iN}\rangle^{2}}$
vs. $L^{-0.5}$ in the DE, GME, and ME distributions along with a linear fit to
$\sigma_{\text{den}}$ versus $L^{-0.5}$. The data strongly suggests that the
fluctuations scale as $L^{-0.5}$ as predicted in Ref. biroli2009  and will
vanish in the thermodynamic limit for all three ensembles.

In the main text, we presented evidence that the fluctuations of the
eigenstate to eigenstate expectation value of the momenta occupations vanish
in the thermodynamic limit for the DE and GME. We have also studied how the
width of the distribution of the momenta eigenstate expectation values,
$\sigma_{k}$, scales for the ME. In Fig. 6(b), we plot the fluctuations of
$\hat{n}_{ka=0}$. For the microcanonical ensemble it is difficult to reach a
conclusion as to whether it vanishes or remains finite in the thermodynamic
limit. We do, however, find that for this observable the results are clearly
different from those for observables that only contain short range
correlations, and that the behavior of the width of its distribution within
the ME behaves quite differently from the one within the DE and the GME. In
the latter two the width decreases much rapidly with system size.

## References

  * (1) M. Rigol,  Phys. Rev. A 72, 063607 (2005).
  * (2) R. Ellis, K. Haven, and B. Turkington, J. Stat. Phys. 101, 999 (2000).
  * (3) M. Costeniuc, R. Ellis, H. Touchette, and B. Turkington, J. Stat. Phys. 119, 1283 (2005).
  * (4) M. Rigol, V. Dunjko, V. Yurovsky and M. Olshanii,  Phys. Rev. Lett. 98, 050405 (2007).
  * (5) G. Biroli, C. Kollath and A. Laeuchli,  ArXiv/0907.3731.


# Discrete Impedance Metasurfaces for 6G Wireless Communications in D-Band

S. Kosulnikov, X. Wang, and S. A. Tretyakov (Corresponding author: X. Wang)S.
Kosulnikov and S. A. Tretyakov are with the Department of Electronics and
Nanoengineering, Aalto University, FI-00079 Aalto, Finland. X. Wang is with
Institute of Nanotechnology, Karlsruhe Institute of Technology, 76131
Karlsruhe, Germany.(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>sergei.tretyakov@aalto.fi).

###### Abstract

Engineering and optimization of wireless propagation channels will be one of
the key elements of future communication technologies. Metasurfaces may offer
a wide spectrum of functionalities for passive and tunable reflecting devices,
overcoming fundamental limits of commonly used conventional phase-gradient
reflectarrays and metasurfaces. In this paper， we develop an efficient way for
the design and implementation of metasurfaces with high-efficiency anomalous
reflector functionalities. The developed numerical method provides accurate,
fast, and simple metasurface designs, taking into account non-local near-field
interactions between array elements. The design method is validated by
manufacturing and experimental testing of highly efficient anomalous
reflectors for the millimeter-wave band.

###### Index Terms:

Metasurface, diffraction grating, reflection coefficient, far-field scattering

##  I Introduction

Novel wireless communication technologies are targeted for the need of growing
data transfer rates, which leads to a shift to higher operational frequencies.
The millimeter-wave (MM-wave) technologies become a basis for new generations
of the wireless communication systems. However, communications at such high
frequencies suffer from high free-space attenuation even in indoor scenarios.
For this reason, it is necessary to use highly directive antennas, losing the
advantages of multi-path propagation that usually provides connectivity even
without a direct line of sight (LOS) between the antennas. One of the very
noticeable techniques potentially solving this issue is based not on
modifications of the source and/or receiver itself, but on modifications and
optimizations of the propagation environment. Indeed, especially indoors, it
is very attractive to realize a scenario where the signal will be efficiently
relayed from some off-site objects. It is worth noting that repeater-like-
devices do not have to be active: even a passive but smart optimization of the
signal propagation path can improve the communication channel by reflecting or
transmitting wave beams towards the desired directions. We also note that this
approach is in line with the “green trends” for efficient and ecological
resource utilization. In this work, we consider a particular example of
passive anomalous reflectors that reflect incident waves to arbitrarily set
reflection angles. This functionality can be used to create an effective
wireless link for non-line-of-sight (NLOS) communications even for MM-wave
systems based on the use of high-gain antennas. Figure 1 illustrates the
considered scenario, where the signal from the transmitter cannot reach the
client receiver in the LOS mode due to a wall obstacle, but still a reliable
wireless channel is set due to redirecting and focusing the reflected signal
by an anomalously reflecting metasurface.

The research of anomalous reflectors has been actively going on from 60’s.
Perhaps the earliest anomalous reflectors are reflectarrays where each element
separated by $\lambda/2$ distance from neighboring elements locally reflects
at different phases, forming a desired reflection wavefront [1, 2]. However,
for grazing reflection angles, the device efficiency significantly drops.
Later, this idea was implanted in arrays with subwavelength unit cells: phase-
gradient metasurfaces [3, 4]. The reflection phase of each unit cell of the
metasurface (realized as a periodical lattice of cells) is linearly and
periodically varying on the surface. Still, such phase-gradient metasurfaces
inherit the main drawbacks of reflectarrays: high reflection efficiency is
achievable only when the deviation from the usual reflection law is not large.
At sharp reflection angles spurious scattering dramatically increases. Since
2016, it is recognized that low efficiency of a device operating far from the
specular or retro-reflection regimes is caused mainly by the impedance
mismatch of the incident and reflected waves [5, 6]. For a given plane-wave
incidence, considering only a single plane wave reflection, it is not enough
to satisfy the local lossless boundary condition of the metasurface [6]. In
[7], it was pointed out that evanescent fields have to be engaged and
optimized to realize a purely reactive boundary that perfectly directs an
incident plane wave to an anomalous direction. Based on this general
principle, different design methods have been proposed. For example, in [8,
9], the evanescent fields are optimized by ensuring local power conservation
on the surface and finding a locally reactive surface. In contrast, the design
methods in [10, 11, 12] are based on the global power conservation where the
structures of meta-atoms are collectively optimized within a super-cell. This
method was recently used for the creation of anomalous reflectors in D-band
[13], where patch sizes in a super-cell were collectively optimized starting
from the local periodic approximation as an initial guess, with good results.
In parallel, the meta-grating method [14, 15] was suggested as another
alternative to realize perfect anomalous reflection with a reduced number of
meta-atoms.

In this paper, we propose an effective design method for the realization of
perfect anomalous reflection into arbitrary directions. As a basis for design,
we select an impedance sheet placed on a grounded dielectric substrate. We
first discretize the impedance sheet into finite elements with uniform values
of the sheet reactance. For this discretized structure we derive an analytical
formula for all the scattered Floquet modes of the metasurface. Then we use an
inverse design optimization method to optimize the discrete impedances until
the amplitudes (and possibly phases) of the scattered modes satisfy the design
objectives. This method is different from the known impedance-based methods
[8, 16], where the sheet impedance is first assumed to be continuous in the
optimization process and discretized into finite elements only after
optimization. After discretization, such metasurfaces may not behave as
perfect as expected from the theoretical results of the continuous-profile
optimization. In contrast, in the proposed method, the sheet impedance is
optimized already in the discretized form. Therefore, the performance does not
suffer from degradations caused by surface discretization. In addition, the
proposed method is not limited by the requirement of the local power
conservation on the surface: the input impedance seen at the input surface can
take purely imaginary or complex values. The overall passivity of the device
is automatically ensured due to the passivity of the impedance sheet and the
substrate material. Therefore, the optimized results encompass both
realization possibilities (local or global power conservation), which is more
general than in the earlier developed methods. In addition, using the same
approach, it is possible to design other devices for wave shaping, such as
multi-channel beam splitters with arbitrarily assigned power ratios.

Figure 1: Application scenario of an anomalous reflector in 6G wireless
communications.

This paper is organized as follows: in the second section we introduce the
theoretical model of the developed numerical optimization method, in the third
section we focus on a particular example implementing an anomalous reflector
operating at D-band, the fourth section is focused on the experimental
validation of the implemented anomalous reflector, the conclusions and
discussion section is closing the paper.

##  II Discrete sheet impedance model for reflection optimization

The metasurface is composed of a discrete impedance array on top of a grounded
dielectric substrate with permittivity $\epsilon_{\rm d}$ and thickness $d$.
Figure 2(a) shows a single unit cell of the metasurface with the spatial
periodicity $D$ (the size of one super-cell). We discrete one super-cell into
$K$ elements, as shown in Fig. 2(a), where the central coordinate of the
$m$-th cell is denoted as $x_{m}={(2m-1)D}/{2K}$. The sheet impedance of the
$m$-th subcell is denoted as $Z_{m}$ where $m\in[1,K]$. The discretized sheet
impedance profile $Z_{\rm s}(x)$ is piece-wise homogeneous and can be viewed
as a cascade of $K$ step functions. As for any periodic function, it can be
decomposed into Fourier series as

$Z_{\rm s}(x)=\sum_{m=-\infty}^{+\infty}g_{m}e^{-jm\beta_{\rm M}x},$ (1)

where $g_{m}$ are the Fourier coefficients. A plane wave illuminates the
metasurface at the angle of $\theta=\theta_{\rm i}$. Due to the periodicity of
the reflector, the scattered field is a sum of an infinite number of Floquet
harmonics, denoted by index $n$. The tangential wavevector of the Floquet
modes scattered from periodical structures can be written as
$k_{xn}=k_{0}\sin\theta_{\rm i}+n\beta_{\rm M}$, where $\beta_{\rm M}=2\pi/D$
is the spatial modulation frequency of the surface impedance, $n$ is the mode
order, $\theta_{\rm i}$ is the incident angle, and
$k_{0}=\omega_{0}\sqrt{\epsilon_{0}\mu_{0}}$ ($\omega_{0}$ is the incident-
wave frequency) is the free-space wavenumber. Modes satisfying
$|k_{xn}|<k_{0}$ propagate into the far-zone at the angles $\theta_{n}$,
defined by

$\sin\theta_{n}=\frac{k_{0}\sin\theta_{\rm i}+n\beta_{\rm M}}{k_{0}}.$ (2)

For $|k_{xn}|>k_{0}$, the modes are evanescent, exponentially decaying along
the surface normal.

Next, based on the mode-matching method, we aim to find the amplitudes of all
scattered modes at the metasurface plane ($z=0$). For convenience of analysis,
we use the transmission-line model shown in Fig. 2(b). The sheet impedance is
modeled as a shunt impedance in the transmission line. The grounded substrate
is modeled as a shorted transmission line with its length equal to the
substrate thickness. The current and voltage in the transmission-line model
are analogous to the total tangential electric and magnetic fields on the
surface ($z=0$). Unlike conventional transmission-line modeling of uniform
metasurfaces, where the current and voltage only have one mode and can be
treated as scalar numbers, in this case, the space-modulated metasurface
excites infinitely many spatial modes (including propagating and evanescent
ones). Therefore, the current and voltage in the transmission line are not
simple scalar numbers but infinite sums:

$I_{\rm s}(x)=\sum_{n=-\infty}^{+\infty}i_{\rm s}^{n}e^{-jk_{xn}x}$ (3a) and
$\quad V_{\rm s}(x)=\sum_{n=-\infty}^{+\infty}v_{\rm s}^{n}e^{-jk_{xn}x}.$
(3b)

Figure 2: (a) Side view of a metasurface containing discrete impedance sheets
on a grounded substrate. (b) Equivalent circuit of the metasurface.

The current and voltage on the impedance sheet must obey Ohm’s law, $V_{\rm
s}(x)=Z_{\rm s}(x)I_{\rm s}(x)$. Substituting Eqs. (3a), (3b), and (1) into
Ohm’s law, we have

$\sum_{n=-\infty}^{+\infty}v_{\rm
s}^{n}e^{-jk_{xn}x}=\sum_{n=-\infty}^{+\infty}\sum_{m=-\infty}^{+\infty}g_{m}i_{\rm
s}^{n}e^{-jk_{x,n+m}x}.$ (4)

By shifting the summation order $n$ to $n-m$ on the right side of Eq. (4), we
get an expression where both sides share the same basis $e^{-jk_{xn}x}$ which
can be eliminated together with the summation over $n$. Therefore, Eq. (4) can
be simplified as

$v_{\rm s}^{n}=\sum_{m=-\infty}^{+\infty}g_{m}i_{\rm s}^{n-m}.$ (5)

One can see from Eq. (5) that every voltage harmonic is related to all current
harmonics due to mode coupling. Considering a finite number of Floquet modes
from $n=-N$ to $n=N$, Eq. (5) can be written in a form of matrix
multiplication:

$\begin{pmatrix}v_{\rm s}^{-N}\\\ v_{\rm s}^{1-N}\\\ \vdots\\\ v_{\rm
s}^{+N}\end{pmatrix}=\begin{pmatrix}g_{0}&g_{-1}&\cdots&g_{-2N}\\\
g_{1}&g_{0}&\cdots&g_{1-2N}\\\ \vdots&\vdots&\ddots&\vdots\\\
g_{2N}&g_{2N-1}&\cdots&g_{0}\end{pmatrix}\begin{pmatrix}i_{\rm s}^{-N}\\\
i_{\rm s}^{1-N}\\\ \vdots\\\ i_{\rm s}^{+N}\end{pmatrix}.$ (6)

Denoting the current and voltage arrays in Eq. (6) as $\mathbf{v}_{\rm s}$ and
$\mathbf{i}_{\rm s}$, and the $(2N+1)$ dimensional impedance matrix as
$\mathbf{Z}_{\rm s}$, Eq. (6) can be written in a simple form:
$\mathbf{v}_{\rm s}=\mathbf{Z}_{\rm s}\cdot\mathbf{i}_{\rm s}$.

Similarly to the sheet impedance, the input impedance of the grounded
substrate can also be written in matrix form. Since the substrate is spatially
uniform and there is no mode coupling in the uniform substrate, the impedance
of the grounded substrate is a $(2N+1)$ dimensional diagonal matrix
$\mathbf{Z}_{\rm gs}$. The $n$-th diagonal term is the input impedance of
grounded substrate for the $n$-th scattering mode (note that we number the
matrix columns and rows from $-N$ to $+N$ instead of from $1$ to $2N+1$):

$\mathbf{Z}_{\rm gs}(n,n)=\frac{1}{Z_{{\rm d},n}\tanh(jk_{zn}^{\rm d}d)}.$ (7)

Here, $Z_{{\rm d}{,n}}$ is the wave impedance of the dielectric substrate
material, and $k_{zn}^{\rm d}=\sqrt{\omega_{0}^{2}\epsilon_{0}\epsilon_{\rm
d}\mu_{0}-k_{xn}^{2}}$ is the normal component of the wavevector in the
dielectric substrate for the $n$-th mode. For TE polarized waves, $Z_{{\rm
d},{n}}=(\mu_{0}\omega_{0})/k_{zn}^{\rm d}$ and for TM polarized waves
$Z_{{\rm d},{n}}=k_{zn}^{\rm d}/(\omega_{0}\epsilon_{0}\epsilon_{\rm d})$.
Details for derivation of (7) can be found in Sec. 4 of the Supplementary
Material of [16].

The total input impedance of the metal-backed metasurface can be calculated as
a parallel connection of the gradient penetrable impedance (characterized by
$\mathbf{Z}_{\rm s}$) of the sheet and the input impedance of the metal-backed
substrate $\mathbf{Z}_{\rm gs}$, i.e., $\mathbf{Z}_{\rm tot}=\mathbf{Z}_{\rm
s}\rm||\mathbf{Z}_{\rm gs}$. After $\mathbf{Z}_{\rm tot}$ is determined, we
can calculate all the scattered modes for a given incidence in terms of a
reflection matrix. The reflection matrix is defined as

$\mathbf{\Gamma}_{\rm TM}=(\mathbf{Z}_{\rm
tot}+\mathbf{Z}_{0})^{-1}\cdot(\mathbf{Z}_{0}-\mathbf{Z}_{\rm tot})$ (8)

for TM-polarized incidence. Here, $\mathbf{Z}_{0}$ is the wave impedance of
free space, which is a diagonal matrix. The $n$-th element of $\mathbf{Z}_{0}$
has the same format with $Z_{{\rm d},n}$, only assuming $\epsilon_{\rm d}=1$.
The reflection matrix for TE-polarized incidence is

$\mathbf{\Gamma}_{\rm TE}=(\mathbf{Y}_{\rm
tot}+\mathbf{Y}_{0})^{-1}\cdot(\mathbf{Y}_{0}-\mathbf{Y}_{\rm tot}),$ (9)

where $\mathbf{Y}_{\rm tot}=\mathbf{Z}_{\rm tot}^{-1}$ and
$\mathbf{Y}_{0}=\mathbf{Z}_{\rm 0}^{-1}$ are the admittance matrices. The
derivation details of Eqs. (8) and (9) are presented in APPENDIX. The
reflection matrix relates the tangential  components of the incident and
reflected fields, i.e., $\mathbf{E}_{\rm r}=\mathbf{\Gamma}_{\rm
TE}\cdot\mathbf{E_{i}}$ for TE waves and $\mathbf{H}_{\rm
r}=\mathbf{\Gamma}_{\rm TM}\cdot\mathbf{H_{i}}$ for TM waves. The incident and
reflected tangential fields ($\mathbf{E}_{\rm i}$, $\mathbf{E}_{\rm r}$,
$\mathbf{H}_{\rm i}$, and $\mathbf{H}_{\rm r}$) should be represented by
$2N+1$-dimensional vertical vectors. The vector contains the complex
amplitudes of the considered harmonics.

The above theory shows how to calculate all the scattered harmonics for a
given set of discrete impedance sheets and illumination waves. However, the
design of metasurfaces is an inverse problem: for a given incidence and
desired reflection harmonics, how to find a proper set of discrete impedance
values? Since there is no analytical solution for this inverse problem, we use
mathematical optimization. Next, we introduce the optimization principles. For
a single TE-polarized plane wave incidence, $\mathbf{E}_{\rm i}$ can be
written as

$\mathbf{E}_{\rm i}=[0,\cdots,0,1,0,\cdots,0]^{T},$ (10)

where the incident mode ($n=0$) is positioned in the middle of the array. For
a known illumination (a given $\mathbf{E}_{\rm i}$), our goal is to find a
proper set of discrete grid impedances of an array that realizes the desired
$\mathbf{E}_{\rm r}$. Here, we focus on the perfect anomalous reflection
functionality with the incident power fully reflected to the $n=+1$ scattering
order. Therefore, the desired reflected field vector can be written as,

$\mathbf{E}_{\rm r}=[0,\cdots,0,0,A_{\rm obj},\cdots,0]^{T},$ (11)

where $A_{\rm obj}=\sqrt{\frac{\cos\theta_{\rm i}}{\cos\theta_{\rm r}}}$ is
the amplitude of reflection that ensures that all the incident power is
directed to the anomalous direction [6]. This goal is realized by optimizing
the discrete impedance values $Z_{1},Z_{2},\cdots Z_{K}$ (purely reactive),
or, equivalently to say, optimizing the reflection matrix
$\mathbf{\Gamma}_{\rm TE}$, until the desired reflection mode amplitude is
maximized in the reflection vector $\mathbf{E}_{\rm r}$. Here, we use the
mathematical optimization tool available in the MATLAB package to find the
optimal values of the impedance sheets. In each trial of the optimization,
MATLAB assumes an array of $Z_{1},Z_{2},\cdots,Z_{K}$ and calculates the
reflected fields. Denoting the calculated amplitude of the desired reflection
mode in each trial as $A_{\rm cal}$, the optimization goal is to find a proper
set of $Z_{1},Z_{2},\cdots,Z_{K}$ that minimizes the cost function defined as

$F(Z_{1},Z_{2},\cdots,Z_{K})=|A_{\rm cal}-A_{\rm obj}|.$ (12)

Employing the MultiStart and $fmincon$ algorithms, MATLAB can search for the
minimum value of $F$ in the multidimensional parameter space.

It is important to mention that this approach is quite general. It can be used
for the design of Floquet metasurfaces with arbitrary power distributions
among all the possible Floquet modes, e.g., beam splitters, only by modifying
the objective function in Eq. (12). In addition, the desired reflection phases
can also be set arbitrarily.

##  III Example of discretized impedance optimization

###  III-A Design goals and optimization results

We target to find designs of anomalous reflecting metasurfaces to realize the
application scenario of Fig. 1 at three potential operational frequencies for
D-band communications: $f_{0}=[144.75;157.75;170.90]$ GHz (the corresponding
operational wavelengths $\lambda_{0}=[2.0725;1.9017;1.7554]$ mm), marked as
Designs 1, 2, and 3 in Table I). The incident TE-polarized plane wave comes
from the normal direction $\theta_{\rm i}=0^{\circ}$, and the goal is to
reflect it into the oblique direction $\theta_{\rm r}=50^{\circ}$. Thus, the
metasurface period is $D=\lambda_{0}/\sin{\theta_{\rm r}}$. We discretize the
impedance sheet into $K=8$ elements, and assume a quartz substrate with the
permittivity $\epsilon_{\rm d}=4.2(1-j0.005)$ and thickness $d=209~{}\mu$m. In
the optimization of the sheet reactances of the elements according to the cost
function Eq. (12), the allowed solutions of $Z_{k}$ were restricted to the
range $[-2000j,+50j]$ Ohm. This constraint is introduced to make the actual
implementations of the array elements easier, because large negative or
positive reactances require some extreme geometries of metal elements that may
be difficult or impossible to fabricate.

As a result of numerical optimization, we get a set of different solutions
minimizing Eq. (12). From them we select the most suitable one with reasonably
decaying amplitudes of the evanescent surface wave harmonics, because high
amplitudes of reactive fields near the array lead to smaller frequency
bandwidth and higher losses. Figure 3 shows an example of optimized mode
amplitudes for the model of an infinite periodic structure. It is shown that
the incident power is fully reflected to the anomalous direction ($n=1$) with
nearly zero specular reflection and $n=-1$ order reflection. The values of the
discretized grid reactances are given in Table I as $Z_{1-8}$ for the
corresponding design frequencies.

TABLE I: Design parameters of the implemented anomalous reflectors. In case of the absence of an element in the sub-cell, the dog-bone element parameter is marked as ”x”. | 1 | 2 | 3
---|---|---|---
$f_{0}$, GHz | 144.75 | 157.75 | 170.90
$Z_{1}$, Ohm | $-319j$ | $-137j$ | $-1229j$
$Z_{2}$, Ohm | $-1686j$ | $-1010j$ | $43j$
$Z_{3}$, Ohm | $-346j$ | $-876j$ | $-1074j$
$Z_{4}$, Ohm | $-138j$ | $43j$ | $-926j$
$Z_{5}$, Ohm | $-991j$ | $-833j$ | $-1250j$
$Z_{6}$, Ohm | $-1721j$ | $-582j$ | $-2000j$
$Z_{7}$, Ohm | $50j$ | $-775j$ | $-141j$
$Z_{8}$, Ohm | $-1140j$ | $-1053j$ | $-2000j$
$C_{x1}$, $\mu$m | 156 | 181.7 | 50
$C_{x2}$, $\mu$m | 64 | 82 | 0
$C_{x3}$, $\mu$m | 205 | 60 | 65
$C_{x4}$, $\mu$m | 286.5 | 0 | 125
$C_{x5}$, $\mu$m | 119.5 | 55 | 155
$C_{x6}$, $\mu$m | 63.7 | 152 | x
$C_{x7}$, $\mu$m | 0 | 150 | 201.7
$C_{x8}$, $\mu$m | 98 | 140 | x
$C_{y1}$, $\mu$m | 40 | 40 | 40
$C_{y2}$, $\mu$m | 40 | 40 | 0
$C_{y3}$, $\mu$m | 40 | 40 | 40
$C_{y4}$, $\mu$m | 40 | 0 | 40
$C_{y5}$, $\mu$m | 40 | 40 | 40
$C_{y6}$, $\mu$m | 40 | 40 | x
$C_{y7}$, $\mu$m | 0 | 40 | 40
$C_{y8}$, $\mu$m | 40 | 40 | x
$L_{x1}$, $\mu$m | 60 | 60 | 60
$L_{x2}$, $\mu$m | 60 | 60 | 116.4
$L_{x3}$, $\mu$m | 60 | 60 | 60
$L_{x4}$, $\mu$m | 60 | 145 | 60
$L_{x5}$, $\mu$m | 60 | 60 | 60
$L_{x6}$, $\mu$m | 60 | 60 | x
$L_{x7}$, $\mu$m | 98 | 60 | 60
$L_{x8}$, $\mu$m | 60 | 60 | x
$L_{y1}$, $\mu$m | 100 | 100 | 40
$L_{y2}$, $\mu$m | 60 | 60 | $\lambda_{0}/10$
$L_{y3}$, $\mu$m | 100 | 60 | 40
$L_{y4}$, $\mu$m | 100 | $\lambda_{0}/10$ | 40
$L_{y5}$, $\mu$m | 60 | 60 | 40
$L_{y6}$, $\mu$m | 60 | 60 | x
$L_{y7}$, $\mu$m | $\lambda_{0}/10$ | 60 | 80
$L_{y8}$, $\mu$m | 60 | 60 | x
$\eta_{\rm eff}^{\rm sim,LL}$, $\%$ | 99.72 | 98.65 | 97.33
$\eta_{\rm eff}^{\rm sim,Lossy}$, $\%$ | 90.78 | 91.81 | 90.76

Figure 3: Amplitudes of scattered harmonics for designed metasurfaces with
optimized step-wise uniform sheets: (a) – Design 1, (b) – Design 2, and (c) –
Design 3.

Finally, we make full-wave simulations with Ansys HFSS to validate the
structure performance on the level of the impedance sheet model, studying a
super-cell with $Z_{\rm g}$ consisting of eight elements at the grounded
dielectric layer. The structure is simulated as an infinite array along $x$
and $y$ with the periodic boundary conditions under TE illumination realized
with a Floquet port excitation. Knowledge of the magnitudes of the Floquet
harmonics for a given periodicity along $x$ found in simulations allows us to
estimate the macroscopic reflection coefficient of the infinite structure [17,
18] and calculate the far-zone field reflected from a finite-size metasurface,
whereas the reflection angle is defined in accordance with the model of
diffraction on infinite periodic structures (2). Typically, the best
numerically optimized solutions can realize the desired anomalous reflection
with the efficiency at the level of 99.99%, if material losses are neglected.

###  III-B Implementation of elements in the super-cell

Our next goal is to find proper geometries of the metallic sub-cells that
would realize the desired grid impedances $Z_{k}$. The grid impedance of a
metallic patch or strip positioned on a grounded substrate can be determined
using the circuit model [19]. To do that, we use the locally periodic
approximation and simulate a single sub-cell using periodic boundary
conditions in both $x$ and $y$ directions along the metasurface plane and find
the reflection coefficient. The grid impedance of the meta-atom can be
determined from the reflection coefficients and the substrate parameters. The
incident angle in the simulation is defined as $\theta_{\rm i}$. The input
impedance of the structure can be expressed in terms of the simulated
reflection coefficient $R$ as

$Z_{\rm in}=\frac{1+R}{1-R}\eta_{0},$ (13)

where $\eta_{\rm 0}=\sqrt{\mu_{0}/\epsilon_{0}}$ is the free-space impedance,
$\mu_{0}$ and $\epsilon_{0}$ are the permeability and permittivity of free
space, respectively [20, 21, 19]. $Z_{\rm in}$ is the input impedance of the
parallel connection of the grid impedance $Z_{\rm g}$ and the grounded
substrate impedance $Z_{\rm gd}$. Next, we extract the grid impedance using

$Z_{\rm g}=\frac{Z_{\rm in}Z_{\rm gd}}{Z_{\rm gd}-Z_{\rm in}},$ (14)

where $Z_{\rm gd}=j\eta_{\rm d}\tan{k_{\rm d}d}$ and $k_{\rm
d}=k_{0}\sqrt{\epsilon_{\rm d}-\sin^{2}\theta_{\rm i}}$ is the propagation
constant in the dielectric substrate with $\eta_{\rm
d}=\eta_{0}/\sqrt{\epsilon_{\rm d}-\sin^{2}\theta_{\rm i}}$ (for TE-polarized
wave) and the free-space wavenumber $k_{0}=\omega\sqrt{\epsilon_{0}\mu_{0}}$
[20].

Figure 4: Geometry of the implemented metasurface: (a) – a sub-cell with
parameters varied in the impedance extraction; (b) – final super-cells
implemented for three targeted designs of perfect anomalous reflectors. (c)
Poynting vector distribution at the $zx$ plane. The shaded circle regions show
the regions of non-local power exchange between free space and the metasurface
volume.

The simplest geometry of capacitive sub-cells is a rectangular patch, but such
elements exhibit very different responses for oblique and normal incidences.
This strong angular dependence leads to some impedance mismatch when the
elements are implemented in the final super-cell, therefore requiring an
additional global optimization of the whole super-cell. For this reason, we
used the so-called dog-bone geometry for the sub-cell metal pattern,
illustrated in Fig. 4(a). This is one example of self-resonant grids [22] that
are known to offer high angular stability of response from high-impedance
surfaces, e.g. [23, 24, 25]. Varying the element’s dimensions we obtain
corresponding parameters of the sub-cell to realize the required impedance
values. Results for the extracted reactance $Z_{\rm g}$ for different
structural parameters for Design 1 are presented in Fig. 5. It is worth noting
a very high level of angular stability comparing the results for the normal
and oblique incidences. The obtained values of the implemented elements’
parameters are presented in Table I, and the top view of the implemented
super-cells is shown in Fig. 4(b). Some of the elements require rather large
negative values of reactance, and our analysis shows that these elements can
be substituted with simply open sub-cells without any element presented. Some
other elements require a small positive reactance which is implemented as
inductive strips with the width $L_{x}$, that is, the length is equal to the
overall $y$ periodicity of the metasurface $L_{y}=\lambda_{0}/10$.

Figure 5: Extracted impedance values and the angular stability analysis.

The final structures are simulated as infinite periodic structures. The
results of simulations in the absence and presence of material losses are
given in Table I as $\eta_{\rm eff}^{\rm sim,LL}$ (with lossless metal) and
$\eta_{\rm eff}^{\rm sim,Lossy}$ (with gold layers of 200 nm thickness),
respectively. The results show that the reflection efficiency slightly ($\sim
1-3$ %) degrades after the implementation of a realistic structure, in
comparison with the performance of the optimized sheet-impedance model. When
metal losses are taken into account, the efficiency degrades towards $\sim 90$
%, mainly due to excited non-propagating surface waves along the metasurface.
The efficiency can be improved if the metal deposition thickness increases
(reducing the Ohmic loss in metal).

Figure 4(c) shows the simulated Poynting vector distribution of Design 1 in
the $xoz$ plane. It is clear that in the regions marked by shaded circles,
there is non-local power exchange between the metasurface volume and free
space. The power flows into the metasurface (virtual loss) and then it is re-
radiated into free space (virtual gain). For the other two implemented designs
we see a qualitatively similar picture showing non-locality of power
reflection. Note that although the given examples shows non-local reflection
and only global power conservation, the developed design method is general and
can find both local and non-local solutions.

##  IV Experimental validation

In order to fabricate the designed anomalous reflectors we used a standard
cleanroom photolithography manufacturing process, with 4-inch quartz wafers
for the substrate. We used AZ514E photoresist, 5 nm Ti and 200 nm Au for both
the front side patterning and the ground plane. This Au layer thickness was
selected in order to reasonably minimize the resistive losses, as they may
significantly affect the device performance.

Figure 6: (a) Schematics and (b) photo of the experiment setup and structure
details of Design 1 under microscope.

For experimental validation, we follow a similar procedure as presented in
[21]. A quasi-optical set-up is used for collimation of the incident
illumination from the horn antenna into a narrow spot on the metasurface
sample under study. This allows us to effectively measure the macroscopic
reflection coefficient as a parameter describing response from the
corresponding infinite structure, avoiding the effects of scattering by the
edges.

The experimental set-up is shown in Fig. 6(a). The used vector network
analyser (VNA) is Keysight PNA Network Analyzer N5225A (10 MHz – 50 GHz) with
WR 5.1 VDI extension units operating at 140 – 220 GHz. A WR-5 rectangular horn
antenna Elmika RHA-015E with 22 dBi directivity was used as the source
antenna. The set-up contains three $90^{\circ}$ parabolic mirrors M1, M2, and
M3 with the corresponding focal distances $F_{1}=25.4$ mm, $F_{2}=152.4$ mm,
and $F_{3}=127$ mm. Note that, the characterization of anomalous reflector
usually needs two-port measurement system, i.e., two horn antennas are
required in the normal and anomalous directions. Here, we use quasi-optical
measurement system and time-gating method to simplify it as one-port system.
The Gaussian beam emitted from the source Horn antenna is collimated by mirror
M1 and converged by M2 towards the normal of sample. The beam is then
reflected to the anomalous direction by the sample. Finally, the beam is
reflected by flat mirror M4 and go back to the horn antenna along the same
route. The efficiency of anomalous reflection is determined by measuring the
reflection coefficient $S_{11}$. Note that in this method, the beam is
reflected twice by the sample. Also, the spurious reflections from all the
mirrors can be corrected by normalization. The normalization is made by
measuring the reflection coefficient when the sample is replaced by a flat
mirror rotated to $\theta_{\rm r}/2=25^{\circ}$. In order to estimate the
retroreflection level for the normal incidence we measured the intensity of
the reflected signal tuning the time-gating parameters for reflection from the
sample. Another anomalously reflected harmonic ($n=-1$) was measured by
upturning the sample. The last measurement is approximate, due to the fact
that the measured area could not be exactly at the same position as the
originally measured for the desired anomalously reflected mode.

Figure 7: Reflection efficiency measured with the quasi-optical method. Curver
marked “$-50^{\circ}$” corresponds to the desired anomalously reflected mode
$n=-1$; Curve “$+50^{\circ}$” is for the same mode measured after flipping the
sample; “$0^{\circ}$ (orient. 1)” curve shows the retroreflection level for
the normal incidence corresponding to the orientation of the sample applied
while measuring the case with“$-50^{\circ}$” and “$0^{\circ}$ (orient. 1)”
corresponds to the same case, but for $+50^{\circ}$.

The experimental results are presented in Fig. 7. Here, the shown data point
gives the efficiency measured at the targeted design frequency. These results
are in good agreement with the theoretically and numerically estimated values
presented in Table I. For all three design cases, one can estimate the
operational band of effective anomalous reflection. As a practically
acceptable level, we define a region where the measured anomalous reflection
efficiency is above 50%. Due to the frequency limit of the used extender (140
GHz) we could not observe the whole operational bands for Designs 1 and 2,
therefore the lower limit of the operational band in these cases is in fact
wider than defined on the experimentally validated data. Design 1 provides
more than half-power anomalous reflection efficiency in the range 140 – 157.4
GHz; Design 2 grants it in the range 140 – 170 GHz; and Design 3 covers the
range 140.6 – 183.1 GHz.

##  V Conclusions

We have presented a method for the design and implementation of metasurfaces
with advanced functionalities based on optimizations of discretized sheet
impedance profiles. The developed fast numerical optimization approach allows
us to directly find useful solutions without the need of further
discretization or detailed numerical tuning of the geometry of the implemented
sub-cells. It is worth noting that the proposed method is general, and not
limited to the considered scenario of anomalous reflection. For example, it
can be used to design perfect anomalous reflectors with larger reflection
angles ($\theta_{\rm r}>50^{\circ}$). By introducing more impedance elements
in one super-cell, it is also possible to realize dynamical beam steering with
high efficiency. Such functionality cannot be realized by meta-gratings with
sparse elements. Moreover, the method can be directly applied to the design of
various beam-shaping surfaces, splitters with phase control of all beams, and
other metasurface devices, just by modifying the optimization objective
function. We consider this method very suitable for the development of
optimized reconfigurable intelligent surfaces because by adding tunable
components to the sub-cells, the beam can be dynamically manipulated.

The considered example of practical implementation was targeted for a
$50^{\circ}$ anomalous reflector operating in D-band, which is one of the
first realizations of anomalous reflectors in this frequency range. The
performance of the designed structure has been first validated with a full-
wave simulation of an infinitely periodic reflector. Performance of the
experimentally realized reflector has been confirmed by quasi-optical
measurements that emulate the response of the infinite structure, allowing
measurements of the macroscopic reflection coefficients for excited Floquet
modes. The experimental results fully agree with the theoretical estimations
and simulations.

##  Appendix A Derivation details

Here, we show derivations of Eqs. (8) and (9). The total tangential electric
and magnetic fields are the sum of tangential incident and reflected fields on
the metasurface plane ($z=0$),

$\mathbf{E}_{\rm tot}=\mathbf{E}_{\rm i}+\mathbf{E}_{\rm r},$ (15)

and

$\mathbf{H}_{\rm tot}=\mathbf{H}_{\rm i}+\mathbf{H}_{\rm r}.$ (16)

Assuming a TE-polarized wave incident on the metasurface, the tangential
magnetic field is related with electric field by the free-space wave
admittance matrix:

$\hat{z}\times\mathbf{H}_{\rm i}=\mathbf{Y}_{0}\cdot\mathbf{E}_{\rm
i},\quad\hat{z}\times\mathbf{H}_{\rm r}=-\mathbf{Y}_{0}\cdot\mathbf{E}_{\rm
r}.$ (17)

Note that the tangential magnetic field flips its sign after reflection,
orienting along $-x$ direction, as shown in Fig. 8(a). The total tangential
fields are related by the impedance boundary condition:

$\mathbf{Y}_{\rm tot}\cdot\mathbf{E}_{\rm tot}=\hat{z}\times\mathbf{H}_{\rm
tot}.$ (18)

Substituting Eqs. (15), (16), and (17) into (18), we can obtain the relation
between the incident and reflected electric fields:

$\mathbf{E}_{\rm r}=(\mathbf{Y}_{\rm
tot}+\mathbf{Y}_{0})^{-1}\cdot(\mathbf{Y}_{0}-\mathbf{Y}_{\rm
tot})\cdot\mathbf{E}_{\rm i}=\mathbf{\Gamma}_{\rm TE}\cdot\mathbf{E}_{\rm i}.$
(19)

Figure 8:  Definitions of the coordinate system and field orientations for TE-
and TM-polarized waves.

For the TM polarization, we define the reflection matrix in terms of
tangential magnetic fields $\mathbf{H}_{\rm r}=\mathbf{\Gamma}_{\rm
TM}\cdot\mathbf{H}_{\rm i}$. In this case, it is more convenient to relate the
tangential electric and magnetic fields via the wave impedance as

$\mathbf{E}_{\rm i}=\mathbf{Z}_{0}\cdot(\hat{z}\times\mathbf{H}_{\rm
i}),\quad\mathbf{E}_{\rm r}=-\mathbf{Z}_{0}\cdot(\hat{z}\times\mathbf{H}_{\rm
r}).$ (20)

Substituting Eqs. (15), (16), and (20) into (18), we obtain

$\mathbf{H}_{\rm r}=(\mathbf{Z}_{\rm
tot}+\mathbf{Z}_{0})^{-1}\cdot(\mathbf{Z}_{0}-\mathbf{Z}_{\rm
tot})\cdot\mathbf{H}_{\rm i}=\mathbf{\Gamma}_{\rm TM}\cdot\mathbf{H}_{\rm i}.$
(21)

## Acknowledgment

This work was supported in part by the European Commission through the H2020
ARIADNE project under grant 871464 and by the Academy of Finland under grant
345178.

We acknowledge the use of facilities of Aalto University at OtaNano –
Micronova Nanofabrication Centre. The authors would like to thank Prof. Sami
Franssila, Dr. Ville Jokinen, Dr. Victor Ovchinnikov, Dr. Mikhail
Omelyanovich, Dr. Andrey Generalov, Dr. Irina Nefedova, and Mr. Pouyan
Rezapoor for their useful comments and discussions related to microfabrication
and mm-wave measurement issues.

## References

  * [1] D. Berry, R. Malech, and W. Kennedy, “The reflectarray antenna,” _IEEE Transactions on Antennas and Propagation_ , vol. 11, no. 6, pp. 645–651, 1963.
  * [2] J. Huang and J. A. Encinar, _Reflectarray Antennas_. Wiley, 2008.
  * [3] D. Pozar, “Wideband reflectarrays using artificial impedance surfaces,” _Electronics Letters_ , vol. 43, pp. 148–149(1), February 2007.
  * [4] 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.
  * [5] N. M. Estakhri and A. Alu, “Wave-front transformation with gradient metasurfaces,” _Physical Review X_ , vol. 6, no. 4, p. 041008, 2016.
  * [6] 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.
  * [7] 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.
  * [8] 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.
  * [9] J. Budhu and A. Grbic, “Perfectly reflecting metasurface reflectarrays: Mutual coupling modeling between unique elements through homogenization,” _IEEE Transactions on Antennas and Propagation_ , vol. 69, no. 1, pp. 122–134, 2020.
  * [10] 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.
  * [11] 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.
  * [12] V. S. Asadchy, A. Wickberg, A. Diaz-Rubio, and M. Wegener, “Eliminating scattering loss in anomalously reflecting optical metasurfaces,” _ACS Photonics_ , vol. 4, no. 5, pp. 1264–1270, 2017.
  * [13] Y. Kato, K. Omori, and A. Sanada, “D-Band Perfect Anomalous Reflectors for 6G Applications,” _IEEE Access_ , vol. 9, pp. 157 512–157 521, 2021.
  * [14] 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.
  * [15] 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\.
  * [16] X. Wang, A. Diaz-Rubio, H. Li, S. A. Tretyakov, and A. Alu, “Theory and design of multifunctional space-time metasurfaces,” _Physical Review Applied_ , vol. 13, no. 4, p. 044040, 2020.
  * [17] A. Díaz-Rubio and S. A. Tretyakov, “Macroscopic modeling of anomalously reflecting metasurfaces: Angular response and far-field scattering,” _IEEE Transactions on Antennas and Propagation_ , vol. 69, no. 10, pp. 6560–6571, 2021.
  * [18] 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.
  * [19] X. Wang, A. Diaz-Rubio, A. Sneck, A. Alastalo, T. Makela, J. Ala-Laurinaho, J.-F. Zheng, A. V. Raisanen, and S. A. Tretyakov, “Systematic Design of Printable Metasurfaces: Validation Through Reverse-Offset Printed Millimeter-Wave Absorbers,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 3, pp. 1340–1351, mar 2018.
  * [20] D. M. Pozar, _Microwave Engineering, 4th Edition_ , 4th ed. Hoboken: John Wiley & Sons, Inc., 2012.
  * [21] X. Wang, A. Díaz-Rubio, V. S. Asadchy, G. Ptitcyn, A. A. Generalov, J. Ala-Laurinaho, and S. A. Tretyakov, “Extreme Asymmetry in Metasurfaces via Evanescent Fields Engineering: Angular-Asymmetric Absorption,” _Physical Review Letters_ , vol. 121, no. 25, p. 256802, dec 2018.
  * [22] I. Anderson, “On the theory of self-resonant grids,” _The Bell System Technical Journal_ , vol. 54, no. 10, pp. 1725–1731, 1975.
  * [23] C. Simovski, P. de Maagt, and I. Melchakova, “High-impedance surfaces having stable resonance with respect to polarization and incidence angle,” _IEEE Transactions on Antennas and Propagation_ , vol. 53, no. 3, pp. 908–914, 2005.
  * [24] A. Vallecchi and A. G. Schuchinsky, “Entwined planar spirals for artificial surfaces,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 9, pp. 994–997, 2010.
  * [25] J. Shabanpour and C. R. Simovski, “Angular and polarization stability of broadband reconfigurable intelligent surfaces of binary type,” _IEEE Access_ , vol. 10, pp. 126 253–126 268, 2022.