1
Supporting Information
for
:
Array
-
Le
vel
I
nverse
D
esign of
B
eam
S
teering
A
ctive
M
etasurfaces
Prachi Thureja
1,
†
, Ghazaleh Kafaie Shirmanesh
1,
†
, Katherine T. Fountaine
2
, Ruzan Sokhoyan
1
,
Meir Grajower
1
, and Harry A. Atwater
1,*
1
Thomas J. Watson
Laboratory of Applied Physics, California Institute of Technology, Pasadena,
California 91125, USA
2
NG Next, Northrop Grumman Corporation, One Space Park, Redondo Beach, CA 90278, USA
†
P. Thureja and G. Kafaie Shirmanesh contributed equally to this work.
*
E
-
mail:
haa@caltech.edu
KEYWORDS
:
inverse design, active metasurface, phased array, beam steering, array
configuration, genetic algorithm,
wavefront engineering
2
1.
Independent scatterer model for
subwavelength antenna arrays
The electric far
-
field calculation, as introduced in equation (1) of the main manuscript, is based
on antenna array theory that was initially developed for phased arrays operating at radio
frequencies (RF).
1
As RF phased arrays can maintain large interantenna spacings without the
introduction of additional diffraction orders, their far
-
field response can be modelled as
that
of a
structured collection of independent scattering elements. However, this assumptio
n becomes non
-
trivial
in the case of
metasurfaces with subwavelength antenna spacings operating at visible or
near
-
infrared wavelengths, as the characteristic period d
x
between neighboring antennas
approaches the electromagnetic near
-
field regime given by
2
d
x
2
/
λ
.
1
Thus, coupling between
neighboring elements can only be neglected for antennas that are spaced at large enough distances
and/or possess strongly confined modes. While the electromagnetic near
-
field of the indium tin
oxide (ITO) based metasurface s
tudied in this work was extensively examined by Kafaie
Shirmanesh
et al.
,
2
we demonstrate here an alternative approach to verify the validity of the
independent scatterer model.
Interantenna coupling causes the actual array phase and amplitude profile to
deviate from the
target design. Consequently, it manifests itself as an increase in the relative magnitude of the
sidelobes due to undesired interference effects. While
near
-
field coupling is
not accounted for in
the analytical model,
it
can be observed in full
-
wave simulations of the entire metasurface
configured with the target array profile. Thus, we perform a comparison of the analytically
calculated radiation patterns to the full
-
wave simulations
2
for forward
-
designed array profiles
(Fig.
S1). We would like to bring the reader’s attention to two apparent effects that are observed in the
cases shown in Fig. S1. Firstly, rather than an increase in the relative magnitude of the sidelobes
(that would indicate interantenna coupling), Fig.
S1a
-
c depict an attenuation of the sidelobes at
broadside angles in the simulated radiation patterns. This effect increases the simulated directivity
D
sim
in comparison to its analytically computed counterpart
D
AF
that is based on array factor
calculations
. Secondly, a strong increase in the relative magnitude of the zero
-
order sidelobe and
thus decrease in
D
sim
is reported for
θ
r
= 70.7° (Fig. S1d). To understand the observed deviations,
we remind ourselves that the analytical calculations were performed f
or omnidirectional scatterers
with
E
antenna
= 1. While this assumption is valid for a broad range of steering angles with the current
nanoantenna design it breaks
at
larger angles. Thus, the reported deviations are attributed to a
simplified antenna model
rather than a manifestation of interantenna coupling.
3
Figure S1. Normalized far
-
field radiation patterns
I/I
max
as a function of the polar angle obtained
through full
-
wave simulations
2
(orange dashed) and analytical array factor calculations (blue). Th
e
results are obtained using forward
-
designed stairstep array profiles for an array of 96 antennas
arranged at a period of
d
x
= 400 nm. The operating wavelength is
λ
= 1510 nm. The simulated
(
D
sim
) and analytically calculated (
D
AF
) directivities at the respective steering angles
θ
r
are: (a)
D
sim
= 51.0 and
D
AF
= 43.5, (b)
D
sim
= 48.3 and
D
AF
= 39.5, (c)
D
sim
= 50.1 and
D
AF
= 39.1, (d)
D
sim
= 22.2,
D
AF
= 23.0.
2.
Phase gradient profiles
Forward
-
designed phase gradient profiles rely on
a constant phase shift
φ
s
between adjacent
antennas. For an incident beam normal to the array,
φ
s
is computed as
3,4
φ
s
= 360°·
d
x
·sin
(
θ
r
)
λ
.
(
S
1)
Here,
d
x
is the characteristic period between neighboring phase antennas,
θ
r
is the steering angle,
and
λ
is the operating wavelength. Wrapping of the phase profiles around 360° allows the design
of blazed grating
-
like structures that steer the reflected beam in t
he desired direction. However,
due to a discrete sampling of the phases at fixed spatial increments
d
x
, the blazed grating of an
antenna array comprises of discontinuous steps, as shown in Fig. S
2
a. The discrete sampling
further results in an aperiodicity
of the phase profiles (Fig. S
2
b), as a complete phase shift of 360°
is not necessarily an integer multiple of
φ
s
. The phase profiles approach periodic blazed gratings
for all steering angles as
d
x
goes to smaller values.
4
Figure S
2
. (a) Phase
φ
as a function of the spatial
x
-
coordinate in an antenna array. The black
dashed line represents an unwrapped, continuous phase profile for steering at
θ
r
= 30°. The blue
curve illustrates the corresponding discretization steps due to sampling at constant
spatial
increments
d
x
= 400 nm. The horizontal grey lines mark the edge values for wrapping around 360°.
The operating wavelength is
λ
= 1550 nm. (b) Phase
φ
vs
. spatial
x
-
coordinate for discretized phase
profile wrapped around 360°. The aperture is
extended to 20 μm to display the aperiodicity of the
wrapped phase profile.
3.
Beam steering performance metrics: directivity
vs.
power efficiency
The figure of merit (FOM) quantifying the beam steering performance in this work was chosen
to be the beam dir
ectivity
D
. It is a unitless quantity that depends on the ratio of the intensity at
the desired steering angle
θ
r
to the amount of power scattered into all directions normalized by the
solid angle, as discussed in equation (2) of the manuscript. Thus, it r
emains unaffected by scaling
of the far
-
field radiation patterns by a constant factor.
Directivity is a common metric used to
analyze the performance of RF phased arrays. An ideal metasurface array with
d
x
= 400 nm
operating at
λ
= 1550 nm (
d
x
/
λ
~ 0.25) approaches performances that are reported with an array of
parallel short dipoles.
5
In addition, the optimized sidelobe level reported in this work corresponds
to values that are generally obtained for phased arrays with a complete phase modulatio
n over
360°.
6
T
he power efficiency
η
is determined by the absolute amount of power that is steered into the
main lobe compared to the total input power. For an array profile with varying amplitudes,
η
is
calculated as
η
(
θ
r
)
=
P
m
(
θ
r
)
P
scat
∙
A
eq
(S2)
where
P
m
is the power scattered into the main lobe steering at
θ
r
and
P
scat
is the total scattered
power. The ratio of
P
m
and
P
scat
is multiplied by the equivalent amplitude
A
eq
that would be
required in an array of antennas with constant amplitude to generate an
equivalent amount of total
scattered power. Thus,
A
eq
=
P
scat
/
P
input
with
P
input
being the input power. Note that
P
input
can be
determined by assuming an ideal reflectarray with constant, unit amplitude and a complete phase
modulation over 360°.
5
Due to th
e strong absorption in the active antenna element, the power efficiency of the beam
steering arrays studied in this work
2
is strongly limited. Consequently, the optimized directivity
case discussed in the manuscript (Fig. S3a
-
b) results a power efficiency
of 0.9%, even though 86%
of the total scattered power is directed into the main lobe. Here, we demonstrate as a proof
-
of
-
concept that the same inverse design algorithm can also be applied to a power efficiency
optimization. For this purpose, the figure of
merit is adapted to
FOM =
η
(
θ
r
).
Figure S3c shows
the optimized array profile as well as the corresponding radiation pattern (Figure S3d) for optimal
power efficiency at
θ
r
= 18.3°. It is to be noted that the increase in power efficiency comes at the
cost
of beam directivity, as the algorithm aims to increase the occurrence of large amplitudes in
the antenna array to enhance efficiency.
Therefore
, the amplitude modulation increases, leading to
a reduction in beam directivity. Meanwhile, the opposite
trend holds true for a directivity
optimization: The inverse design aims to minimize amplitude modulation to reduce sidelobes. As
the main phase shift occurs in a low amplitude regime, the minimization of amplitude modulation
results in reduced power effic
iencies. For reference, the corresponding directivity and efficiency
values are tabulated in Table S1.
Figure S3. Inverse
-
designed array amplitude (black dotted) and phase (blue) profiles and far
-
field
radiation patterns for an optimization of
directivity (a)
-
(b) and an optimization of power efficiency
(c)
-
(d), respectively. Optimizations are performed for the ITO
-
based metasurface studied in this
work.
2
The corresponding directivity and power efficiency values are listed in Table S1.
Directiv
ity
D
Efficiency
η
Forward design, stairstep
39.5
2.1%
Inverse design, directivity opt.
72.7
0.9%
Inverse design, efficiency opt.
41.9
2.7%
Table S1. Directivity
D
and efficiency
η
, respectively, for the electro
-
optically tunable, ITO
-
based
metasurface introduced by Kafaie Shirmanesh
et al.
2
The corresponding values are compared for
three different cases: forward
-
designed stairstep array profiles, and inverse
-
designed profiles
optim
ized for directivity and efficiency, respectively.
6
As the scattered light amplitudes are the limiting factor for power efficiencies in beam
steering metasurfaces, we would like to remark that they can be strongly enhanced with the use of
active metasurfa
ces exhibiting higher reflectance / transmittance values, such as all
-
dielectric
metasurfaces.
7
-
9
4.
Complex
dielectric
permittivity of indium tin oxide (ITO)
The continuity of the normal electric displacement component at the interface between two
media
requires
ε
1
∙
E
1
⊥
=
ε
2
∙
E
2
⊥
(
S3
)
where
ε
i
is the complex dielectric permittivity and
E
i
⟂
퐸
푖
⊥
is the normal electric field component
in medium
i
. Hence, operation at an epsilon
-
near
-
zero (ENZ) condition results in a strong field
enhancement in the active ITO layer. The spectral overlap of this ENZ transition with the magnetic
dipole resonance of the antenna thus ensures a strong modulation of the
scattered light respons
e.
The spatial variation of the real and imaginary parts of the ITO permittivity (
ε
ITO
)
under
applied bias are presented in Figure S
4
.
The ITO properties are chosen as described in Ref. 10.
As
can be seen in Fig. S
4
a, for a sufficiently large applied bias, an ENZ condition holds in the ITO
layer where the real part of the ITO permittivity can take values between
-
1 and +1. Figure S
4
b
sho
ws that the imaginary part of the ITO permittivity in the mentioned regions takes nonzero
values. The nonzero complex permittivity of the ITO layer in the ENZ region leads to a finite
electric field confinement in the accumulation region of the ITO. Nevert
heless, owing to the near
-
zero real part of the permittivity, the ITO layer experiences a large field enhancement, as shown
in the electromagnetic near
-
field distributions provided by Kafaie Shirmanesh
et al
.
2
in Parts 2 and
3 of the Supporting Information
of their manuscript
.
Figure S4. Spatial variation of complex dielectric permittivity of ITO under applied bias. (a) Real
and (b) imaginary part of the permittivity of ITO as a function of position for different applied bias
voltages at the operating wavelength of
= 1510 nm.
The gray
-
shaded region in (a) shows the
ENZ regime where
-
1 < Re{
ε
ITO
} < +1.
z
= 5 nm denotes the interface of the ITO and the HAOL
gate dielectric, and
z
= 0 nm represents the interface of the ITO and the Al
2
O
3
dielectric layer
, as
indicated in the inset of (b)
.
Voltage is applied between the ITO and the top Au layer.
7
Figure S5. Optical response of the metasurface incorporating an artificial ITO layer
with zero
collision frequency, and hence zero imaginary part of permittivity. (a) Spatial variation of the real
part of the permittivity of the artificial ITO film as a function of position for different applied bias
voltages at the wavelength of
= 1510
nm. Here,
z
= 5 nm denotes the interface of the ITO and
the HAOL gate dielectric, and
z
= 0 nm represents the interface of the ITO and the Al
2
O
3
dielectric
layer. (b) Reflectance and (c) phase of reflection spectra for different applied biases. (d)
Reflectance and phase of the reflection as a function of applied bias at the operating wavelength
of
= 1490 nm.
According to the Drude model, the co
mplex permittivity of the ITO can be formulated as
ε
ITO
(
ω
)
=
ε
∞
-
ω
p
2
ω
2
+
i
ω
Γ
(
S4
)
where
ε
∞
is the infinite frequency permittivity,
ω
is the angular frequency,
ω
p
is the plasma
frequency and
Γ
is the collision frequency. The latter contributes to Im(
ε
ITO
) and quantifies the
losses in form of absorption through electron
-
electron collisions. Thus, for increased values of
Im(
ε
ITO
), lower scattered light amplitudes are expected as a result of enhan
ced absorption.
1
1
Here,
we simulate the same metasurface structure with an artificial ITO layer with zero imaginary part.
To this end, we set the collision frequency of the ITO layer to be zero (
Γ = 0
in eq. (
S4
)). Figure
S
5
a presents the real part of the
permittivity of the artificial ITO layer as a function of position
(with
z
= 5 nm being the interface of the ITO and the HAOL gate dielectric) for different applied
bias voltages. For sufficiently large applied bias voltages, an ENZ condition holds in the
artificial
8
ITO layer. The spectra of the reflection amplitude (reflectance) and phase are depicted in Fig.
5
b
and c, respectively. Upon changing the applied bias, a notable reflectance and phase modulation
could be obtained. Figure S
5
d illustrates the refl
ectance and phase values as a function of applied
bias at the operating wavelength of
= 1490 nm. As can be seen, a remarkable amplitude
modulation with increased amplitude values at negative voltages, accompanied by a phase shift of
320
o
can be achieved
via
incorporating the mentioned artificial ITO layer within our metasurface.
Notably, reflectance is reduced at ENZ condition even though the active layer is lossless. This
implies that Au absorbs more light due to a stronger field confinement in the gap.
Furthermore,
the resonant position has been shifted
due to the zero imaginary part and a slight change in the real
part of the permittivity of the artificial ITO layer compared to the original ITO film.
Finally, w
e
would like to note that while there have
been theoretical studies on active metasurfaces that are
based on low
-
loss materials
with near
-
zero complex dielectric permittivity
, such as cadmium oxide
(CdO),
1
2
they have yet to be demonstrated experimentally
.
5.
Forward designs in non
-
ideal
antenna arrays
Device non
-
idealities for active metasurfaces include tunable antennas with (i) non
-
unity
amplitudes, (ii) reduced phase modulation range, and (iii) covarying phase and amplitude values.
While non
-
unity amplitudes impact the power efficienci
es of steered beams, the latter two device
characteristics are directly translated into limitations of forward designs, as discussed below. The
phase and amplitude profiles shown in Fig. S
6
are obtained with the scattered light properties of
the electro
-
op
tically tunable metasurface in Ref
.
2
.
Reduced phase modulation range
A reduced phase modulation range in an active metasurface requires modification to the
phase gradient profile to ensure that it remains within the maximal acquired phase shift.
Here, we
discuss two possible adjustments to ideal forward designs. Linear truncated phase profiles consist
of blazed gratings that are truncated symmetrically around 180°. Phase profiles are then shifted
such that the minimal acquired phase value is 0° (F
ig. S
6
a, blue). Step profiles, on the other hand,
simplify the design of phase gradient profiles by repeating a discrete number of phase values over
several antenna (Fig. S
6
c, blue). As the truncated linear phase profiles have a closer resemblance
to ideal
blazed gratings, higher directivities are reported in that case.
Covariation of amplitude and phase
Covarying amplitude and phase values inhibit the design of pure phase gratings, as
indicated by the black dotted lines in Fig. S
6
a, c. As a consequence, i
ncreased destructive
interference results in additional scattering that appears in form of undesired sidelobes. Figure S
6
b,
d illustrate how the far
-
field radiation patterns change after consideration of the phase
-
amplitude
interdependencies for steering a
t
θ
r
= 18.3°. For the two cases considered here, directivities drop
from
D
lin,const
= 72.1 to
D
lin,covary
= 50.7 for linear phase profiles, and from
D
step,const
= 54.8 to
D
step,covary
= 39.5 for stairstep phase profiles.
9
Figure S
6
. (a) Linear truncated p
hase (blue) and corresponding amplitude profile (black dotted)
over 96 antennas for steering at
θ
r
= 18.3°. (b) Normalized far
-
field intensity
I/I
max
vs
. polar angle
θ
for linear truncated phase profiles with constant amplitudes (violet dashed) and covaryi
ng
amplitudes (orange). (c) Stairstep phase (blue) and corresponding amplitude profile (black dotted)
over 96 antennas for steering at
θ
r
= 18.3°. The stairstep phase profile is obtained by periodically
repeating each phase shift of 270°, 180°, 90°, 0° ove
r 3 consecutive antennas. (d) Normalized far
-
field intensity
I/I
max
vs
. polar angle
θ
for stairstep phase profiles with constant amplitudes (violet
dashed) and covarying amplitudes (orange). The phase
-
amplitude relation is obtained from the
optical response of the electro
-
optically tunable metasurface i
ntroduced in Ref. 2
. Antennas are
arr
anged at a period of
d
x
= 400 nm. The operation wavelength is
λ
= 1510 nm. Radiation patterns
are computed with the assumption of omnidirectional antennas.
6.
Iterative genetic optimization
: numerical framework
Figure S
7
outlines the implementation of iterati
ve genetic algorithms using the global
optimization toolbox on MATLAB. The input of the algorithm comprises of the steering angle
θ
r
,
as well as the objective function
FOM
(
x
,
φ
(
V
),
A
(
V
)) that takes into account the tunable scattered
light properties of the
metasurface.
x
is the 1D vector representing the array configuration that
needs to be optimized. In addition, we define the following global variables: the total number of
antennas
N
tot
, the number of optimization rounds
r
tot
, as well as an array containi
ng the number of
possible variables
nvars
which are to be optimized in each iteration. For the active metasurface
with 96 tunable antennas,
nvars
is defined as [4, 8, 24, 48, 96] such that the optimal solution is
found within a maximal number of five itera
tions. The concept of iterative genetic optimization
relies on an initially reduced search space.
10
Figure S
7
. Flowchart of iterative genetic optimization for an array of
N
tot
= 96 antennas. The inner
loop represents the iterative genetic optimization with increasing variable size to approach the high
dimensionality of the underlying problem. The outer loop describes a series of optimization rounds
that allow to take the optima
l solution over multiple repetitions. The latter is required due to the
stochastic nature of genetic optimization.
The algorithm aims to optimize for a sequence of small number of variables that are
periodically repeated over the entire array. Once an optimal solution
x
opt
is found,
nvars
is
incrementally increased to the next value. An initialization with
k
=
nvars
(
i
+
1)/
nvars
(
i
) repetitions
of the current optimized solution
x
opt
guides the algorithm in larger solution domains. This
procedure is repeated until all
N
tot
antennas are considered as free variables in the final iteration.
Once
nvars
(
i
) =
N
tot
, the current op
timization round is terminated and
x
opt
is stored along with its
corresponding function value in an array. This iterative optimization process is repeated for
r
tot
rounds, after which the solution with the maximal
f
opt
is given as output. This step is nece
ssary due
to the stochastic nature of genetic optimization. Note that prior knowledge from blazed grating
design allows us to make the algorithm more efficient. The number of variables that are to be
optimized in the first iteration can be determined as a
function of the steering angle
θ
r
, using the
grating equation defined in (1).
11
Convergence properties
The iterative optimization process for the active metasurface analyzed in our study is
illustrated in the supporting movie ‘MovieS1.avi’.
The target steering angle is
θ
r
= 18.3°.
As can
be seen in the convergence plot, t
he main contribution to the increase in directivity occurs in the
initial optimization in the reduced solution space. Once an optimal solution is found, the algorithm
moves
on to the next bigger s
olution domain
where minor changes in phase and amplitude lead to
a reduction of the sidelobes
which
are displayed
in the log
-
scaled inset of the
far
-
field radiation
pattern
as well as in the best directivity value
.
Since a separate
optimization was performed in
order to collect the required data at each generation, the optimized directivity differs from the value
reported in the manuscript
.
Computational cost
In contrast to forward
-
designed array profiles that rely on an analytical equation, as discussed
in (S1), the inverse design approach
comes with
enhanced computational cost due to a
consideration of the antenna
-
specific functional response. For the problem
s analyzed in our work,
t
he optimal solution
in
each iteration is generally obtained within 200
-
600 generations. A single
computation (
r
tot
= 1) for an iterative optimization of 96 variables performed on our workgroup
computer (
Intel Xeon E5
-
2687W processo
r, 20 cores
) takes approximately 12 min.
The required computation time highly depends on the total number of variables that are to be
optimized. Figure S
8
shows the average computation time TOC
avg
for a single iterative
optimization round as a function of
N
tot
. The average time was evaluated over
r
tot
= 10 optimization
rounds for six different steering angles (
θ
r
= 9.0°, 10.9°, 13.6°, 18.3°, 28.1°, 70.7°). Notably, the
computation time scales linea
rly as
O
(
N
tot
) in the investigated regimes while the solution space
scales exponentially as
O
(
s
N
tot
) where
s
is the number of sampling points. For our study, the
antenna
-
specific
scattered light response was sampled at
s
= 65 discrete voltage points. The
d
ifference in scaling is attributed to the stopping criteria:
In
the current implementation
, the
algorithm stops once the average change in best function value over 250 generations is less than
10
-
6
. As the most significant contribution of the directivity e
nhancement occurs for the initial
optimization in a reduced solution domain, each subsequent iteration adds approximately 250
generations to the optimization
process
that result in minor performance enhancements
. Therefore,
a linear increase in computation
time is observed.
This phenomenon can
also
be seen in the
supporting movie ‘MovieS1.avi’, illustrating the convergence of the
iterative optimization.
In
future work, the stopping criteria can be optimized such that the computational cost is reduced
withou
t a significant loss in best performance.