# Optimal Inverse Design Based on Memetic Algorithms – Part 2: Examples and
Properties

Miloslav Capek, , Lukas Jelinek, Petr Kadlec, and Mats Gustafsson Manuscript
received ; revised . This work was supported by the Czech Science Foundation
under project No. 21-19025M.M. Capek and L. Jelinek are with the Czech
Technical University in Prague, Prague, Czech Republic (e-mails:
{miloslav.capek; lukas.jelinek}@fel.cvut.cz).M. Gustafsson is with Lund
University, Lund, Sweden (e-mail: mats.gustafsson@eit.lth.se).P. Kadlec is
with the Brno University of Technology, Brno, Czech Republic (e-mail:
kadlecp@feec.vutbr.cz).Color versions of one or more of the figures in this
paper are available online at http://ieeexplore.ieee.org.

###### Abstract

Optimal inverse design, including topology optimization and evaluation of
fundamental bounds on performance, which was introduced in Part 1, is applied
to various antenna design problems. A memetic scheme for topology optimization
combines local and global techniques to accelerate convergence and maintain
robustness. Method-of-moments matrices are used to evaluate objective
functions and allow to determine fundamental bounds on performance. By
applying the Shermann-Morrison-Woodbury identity, the repetitively performed
structural update is inversion-free yet full-wave. The technique can easily be
combined with additional features often required in practice, e.g., only a
part of the structure is controllable, or evaluation of an objective function
is required in a subdomain only. The memetic framework supports multi-
frequency and multi-port optimization and offers many other advantages, such
as an actual shape being known at every moment of the optimization. The
performance of the method is assessed, including its convergence and
computational cost.

###### Index Terms:

Antennas, numerical methods, optimization methods, shape sensitivity analysis,
structural topology design, inverse design.

##  I Introduction

Inverse design is a time-consuming process with no certainty regarding global
minimum feasibility [1]. This is generally the case, no matter how
sophisticated the method employed [2, 3, 4, 5, 6, 7, 8, 9] and there is no
proof of convergence towards the global minimum [2, 3]. However, the global
minimum is typically not needed in practice. A sufficiently good solution has
to be found in a reasonable time. As such, a good balance between detailed
local search and large-scale exploration of the solution space has to be
achieved [10]. These properties are provided by the approach introduced in
Part 1 [11] which lays down the groundwork for this second part.

The inverse design procedure proposed in this paper relies on two equally
important steps. In the first step, the fundamental bound on an optimized
metric is evaluated. This provides a stopping criterion for the memetic scheme
and delimits the performance of any realized device. The memetic scheme is
initiated in the second step and attempts to minimize the distance from the
fundamental bound. It may happen in some cases that the bound is not
reachable, e.g., there are not enough feeders, or they are not properly
placed. In such a case, the optimization can be restarted with improved
initial conditions.

The memetic optimization method combines two distinct approaches – local and
global steps. The local step is based on investigating the smallest topology
perturbations, i.e., changes in the value of an objective function if one of
the degrees-of-freedom (DOF) is removed or added to the optimized structure.
The local step was proposed in [12] and [13], where only DOF removals were
used to detect the local minima. The approach was extended by the possibility
of adding DOF to the system in [14]. Satisfactory performance of the local
step was confirmed on Q-factor minimization [14], as well as on the
minimization of reflectance of a pixel antenna [15] or, recently, the design
of surface unit cell [16]. Part 1 [11] merged both smallest perturbations
(addition and removal of DOF) into a unified framework and combined it with
the global step.

Since a fixed discretization grid is used, the differences calculated from the
change of the objective function value under the smallest topology
perturbations (topology sensitivities) represent a discrete analogue to
gradient over structural variables. As with gradient-based optimization
schemes, these differences are used to search for a local minimum via an
iterative greedy search [17].

The global step is designed to restore and maintain diversity when the local
step finds a local minima. Heuristic approaches are known for their robustness
[18, 19], applied, e.g., in the form of genetic algorithms operating over
locally optimal shapes. The binary nature of genetic algorithms suits the
combinatorial-type optimization solved in this work. While heuristics do not,
in general, perform well [20], only good properties, including versatility,
robustness, and easy implementation, are used here while the disadvantages
(mainly computational requirements and slow convergence [21]) are mitigated by
the underlying local step. Moreover, it is shown in this paper that, in many
cases, only the local step is needed to identify shapes good enough for
practical purposes. The above-mentioned properties and claims are confirmed in
this paper using four examples involving electrically small and medium-sized
problems. Both scattering and antenna scenarios are treated. The minimizing
functions are single- and multi-objective.

The paper is structured as follows. Implementation and benchmarking details
are provided in Section II. Section III deals with the minimization of the
Q-factor. Results for a rectangular plate and a spherical shell are
consistently compared to fundamental bounds. The influence of discretization
on precision and computational time is provided. Section IV generalizes the
objective function to a multi-objective case studying the trade-off between
Q-factor and input impedance (matching). A code implementing the local step
for minimizing Q-factor is published as supplementary material [22]. An array
is synthesized in Section V to maximize realized gain. The different number of
elements and different spacings are considered. Section VI describes how to
optimize a region adjacent to a lossy chip in which the power absorbed from
the incoming plane wave has to be maximized. Finally, the various aspects of
the method are thoroughly discussed in Section VII and the paper is concluded
in Section VIII.

##  II Examples – Methodology

The properties and overall performance of the proposed optimization procedure
are discussed using various examples focusing on different aspects of the
method. All parts were implemented in MATLAB [23], the matrix operators were
evaluated in AToM [24], and the genetic algorithm from FOPS [25] was utilized.
The surface method of moments (MoM) for good conductors was utilized with Rao-
Wilton-Glisson (RWG) basis functions [26] defined over a Delaunay
triangulation [27]. The $5$-th order quadrature rule [28] was applied to
evaluate MoM-based integrals.

The examples with recorded computational time were evaluated on a computer
with an AMD Ryzen Threadripper 1950X CPU (16 physical cores, 3.4 GHz) with 128
GB RAM. The extensive studies, sweeping one or more parameters, were evaluated
on an RCI cluster [29]. All parts of the code were implemented as described in
[11] with the exception of the GPU evaluation of an objective function, i.e.,
the local step was vectorized. The global step is calculated with the help of
parallel computing.

##  III Minimum Q-factor

Minimization of Q-factor has a long history [30, 31] and is still a topical
problem. The importance of this quantity for electrically small antennas is
given by its inverse proportionality to fractional bandwidth [32], a parameter
suffering greatly from small electrical size [33]. While the determination of
the lower bound to quality factor has matured [34], the corresponding optimal
shapes are known thanks only to empirical designs [35, 36, 31]. This
inevitably limits their exclusive applicability to cases when they are known
to be approximately optimal (e.g., a folded spherical helix for a spherical
region [37] or a meanderline for a rectangular region [31]).

Here, we demonstrate how to utilize the technique presented in Part 1 of this
paper [11] to get solutions close to the fundamental bounds in an automated
manner. The optimization setup is chosen so that the results can be compared
with empirically synthesized structures [31], i.e., with the region
$\varOmega$ in the form of a rectangle as, e.g., described in Fig. 1 of Part 1
[11].

###  III-A Inverse Design

The rectangular region of varying mesh grid density and varying electrical
size expressed in $ka$, where $k$ is the wavenumber and $a$ is the radius of
the smallest circumscribing sphere, is considered first. The optimization was
performed with $N_{\mathrm{A}}=48$ agents (three times the number of physical
cores of the Threadripper 1950X processor) utilized during the global step
[11]. The excitation is realized via a discrete delta gap feeder placed in the
middle of the longer side, close to the boundary of the design region. The
amplitude and phase of the excitation does not affect the result thanks to the
quadratic forms used to define the Q-factor as

$\begin{split}Q(\mathbf{g})&=Q_{\mathrm{U}}(\mathbf{g})+Q_{\mathrm{E}}(\mathbf{g})\\\
&=\dfrac{1}{2}\dfrac{\mathbf{I}^{\mathrm{H}}\left(\mathbf{g}\right)\mathbf{W}\mathbf{I}\left(\mathbf{g}\right)}{\mathbf{I}^{\mathrm{H}}\left(\mathbf{g}\right)\mathbf{R}_{0}\mathbf{I}\left(\mathbf{g}\right)}+\dfrac{1}{2}\dfrac{\left|\mathbf{I}^{\mathrm{H}}\left(\mathbf{g}\right)\mathbf{X}_{0}\mathbf{I}\left(\mathbf{g}\right)\right|}{\mathbf{I}^{\mathrm{H}}\left(\mathbf{g}\right)\mathbf{R}_{0}\mathbf{I}\left(\mathbf{g}\right)}\end{split}$
(1)

here decomposed into its untuned $Q_{\mathrm{U}}$ and tuning $Q_{\mathrm{E}}$
parts [11]. This corresponds to the classic Q-factor definition [32]

$Q=\dfrac{2\omega\max\left\\{W_{\mathrm{m}},W_{\mathrm{e}}\right\\}}{P_{\mathrm{rad}}}=\dfrac{\omega
W}{P_{\mathrm{rad}}}+\dfrac{\left|P_{\mathrm{react}}\right|}{2P_{\mathrm{rad}}},$
(2)

where $W,W_{\mathrm{m}},W_{\mathrm{e}}$ are cycle mean values of the total,
magnetic and electric stored energy, respectively, and where
$P_{\mathrm{react}},P_{\mathrm{rad}}$ are the reactive and radiated powers,
see Appendices B and C of Part 1, [11]. The materials allowed by the optimizer
are a perfect electric conductor (PEC) and vacuum.

Figure 1: Shape optimization (inverse design) of an antenna minimizing
Q-factor. The rectangular plate is used as the region $\varOmega$, see
examples for grids $8\times 5$ and $16\times 10$ as the insets in the top-
right corner. Due to cubic dependence on electrical size, Q-factor values are
normalized to $(ka)^{3}$, [31]. The optimized performance for four mesh grids
(different line styles) is compared with fundamental bounds $Q_{\mathrm{lb}}$
and $Q_{\mathrm{lb}}^{\mathrm{TM}}$, [34, 38]. The self-resonance test
$Q_{\mathrm{E}}/Q_{\mathrm{U}}$ is depicted in the top pane. The number of
meanders which are required to construct self-resonant meanderline antenna
reaching the bound $Q_{\mathrm{lb}}^{\mathrm{TM}}$ and fitting the design
region is indicated by diamond marks at the top. The number of meanders is
inferred from the parameterization used in [31].

The algorithm was set as follows: the local step can potentially have
infinitely many iterations $I$, however, it is terminated when the relative
difference $\epsilon_{\mathrm{loc}}=10^{-7}$ between two consecutive
iterations is reached. Similarly, the relative difference for the global step
was set to $\epsilon_{\mathrm{glob}}=10^{-7}$ with a maximum of $J=250$ global
iterations.

###  III-B Comparison with the Fundamental Bounds

The results are compared with bounds $Q_{\mathrm{lb}}$ and
$Q_{\mathrm{lb}}^{\mathrm{TM}}$, shown in Fig. 1, and confirm the good
performance of the method. The placement of the feed and reflection symmetry
of the optimized region makes it possible to excite purely TM-modes [38],
therefore, $Q_{\mathrm{lb}}^{\mathrm{TM}}$ is a tighter bound as compared to
the $Q_{\mathrm{lb}}$ bound. The $Q_{\mathrm{lb}}^{\mathrm{TM}}$ bound is
almost reached in the $ka\in[0.6,\,1]$ region. The top pane of Fig. 1
demonstrates that the structures meeting the bound are close to self-
resonance, i.e., $Q_{\mathrm{E}}/Q_{\mathrm{U}}\rightarrow 0$, following the
expectation that Q-factor reaches its minimum for a self-resonant current
[39]. However, there are two regions where the bound is not met, namely,
$ka<0.6$ and $ka>1$.

The region $ka<0.6$ requires finer granularity of the optimized domain than
what was used in this example. This can be achieved only at the cost of
increasing the number of unknowns. This is verified by the study in [31],
where the number of meanders required to construct self-resonant Q-factor-
optimal meanderlines approximately increases with $(ka)^{-2}$, see [31, Fig.
3]. The necessary number of meanders to reach the bound is shown by orange
diamond markers in Fig. 1. It can be verified that the bound is reachable for
a mesh grid of $20\times 12$ pixels at $ka=0.6$ where approximately $9$
meanders would be needed. This is the maximum number of meanders that this
mesh grid can accommodate.

The region $ka>1$ is meshed densely enough, as can be seen from a comparison
of the black lines in this region. However, the problem is the fixed number of
feeders and their placement. It is highly probable that more than one feeder
and asymmetric placements are needed to reach the lower bound on Q-factor111It
should be noted that (for given material support) to every optimal current
vector $\mathbf{I}_{\mathrm{opt}}$ there exists an optimal vector of
excitations $\mathbf{V}_{\mathrm{opt}}=\mathbf{Z}\mathbf{I}_{\mathrm{opt}}$. A
study provided in [40] further shows that even for a lower number of
controllable excitations, an optimally excited current density on realized
structure attempts to mimic that of the fundamental bound..

Figure 2: Evolution of optimized shapes with global iteration $j$ for a mesh
grid of $12\times 7$ square pixels. Two electrical sizes, $ka=0.66$ and
$ka=1.03$ are treated cf., Fig. 1. Only $N_{\mathrm{A}}=16$ agents were used
for the global step and each picture depicts their average. A level of opacity
(a value between zero and one) equals the average gene over the pool of agents
of the given generation $j$, i.e., the more opaque the color is, the more
agents have this DOF enabled. This postprocessing technique is used only to
show the diversity of the agents for a given global iteration $j$. The last
column corresponds to the final global iteration ($J=38$ for $ka=0.66$ and
$J=35$ for $ka=1.03$). The solid red lines represent basis functions which are
enabled, i.e., current paths enabled by the memetic algorithm.

Two optimized shapes for $ka=0.66$ with $(ka)^{3}Q\approx 5.44$ and
$Q_{\mathrm{E}}/Q_{\mathrm{U}}\approx 0.12$, and $ka=1.03$ with
$(ka)^{3}Q\approx 5.52$ and $Q_{\mathrm{E}}/Q_{\mathrm{U}}\approx 2.3\cdot
10^{-6}$ are depicted in Fig. 2, right column. The second and the fifth global
iterations are shown as well, see the left and middle columns. In all cases,
the presence of a given DOF was evaluated as the average of all
$N_{\mathrm{A}}=48$ agents at a given $j$-th global iteration. No post-
processing of the shapes was performed to show the raw results of the
optimization algorithm. However, additional regularity constraints might be
imposed to improve the manufacturability of the results [41], e.g., the
topology sensitivity map of the final design can be used to manually remove
parts of the structure that obstruct manufacturing but have only minuscule
impact on the optimized metric. It is seen that, while all agents at all
iterations are always local minima (due to the performance of the local step),
there is still high uncertainty about the final structure at the beginning,
while the global minimum becomes uniquely determined close to the end of the
optimization. Notice also, that high uncertainty at the beginning (each agent
sits in different local minima) means high diversity, which is required for
non-convex problems. On the other hand, the technique converges quickly, as
shown in the next section.

###  III-C Computation Time

The performance of the optimization technique has also to be judged in terms
of computation time and convergence rate. The optimization from Fig. 1 is
repeated for $ka=0.5$ for varying mesh grids, see Table I. Depending on the
number of DOF $N$, the number of required global iterations $J$ is shown,
together with the total required computational time $t(J)$, distance from the
bound $Q/Q_{\mathrm{lb}}^{\mathrm{TM}}$, and number of investigated shapes for
removals ($\sum_{i}|{\mathcal{R}}(\mathbf{g}_{i})|$) and additions
($\sum_{i}|{\mathcal{A}}(\mathbf{g}_{i})|$). The last two columns represent
the number of antennas evaluated with the full-wave exact reanalysis method
during the optimization. The number of optimized variables is $B=N-1$ since
one DOF is used for (fixed) delta gap feeding. It is seen that it reaches up
to hundreds of millions of antenna shapes evaluated with (full-wave) MoM.

TABLE I: Total number of global iterations $J$ used, computational time $t$, reached Q-factor, and number of shapes evaluated for removals and additions based on the number of DOF $N$. The values were recorded for settings from Fig. 1 and electrical size $ka=0.5$ grid | $N$ | $J$ |  $t$ (s) | $\dfrac{Q}{Q_{\mathrm{lb}}^{\mathrm{TM}}}$ | $\sum\limits_{i}|{\mathcal{R}}(\mathbf{g}_{i})|$ | $\sum\limits_{i}|{\mathcal{A}}(\mathbf{g}_{i})|$
---|---|---|---|---|---|---
$8\times 5$ | $227$ | $27$ | $24$ | $1.16$ | $8.75\cdot 10^{5}$ | $8.78\cdot 10^{5}$
$12\times 7$ | $485$ | $49$ | $213$ | $1.13$ | $3.86\cdot 10^{6}$ | $4.02\cdot 10^{6}$
$16\times 10$ | $934$ | $85$ | $2190$ | $1.07$ | $1.64\cdot 10^{7}$ | $1.84\cdot 10^{7}$
$20\times 12$ | $1408$ | $146$ | $10400$ | $1.06$ | $3.85\cdot 10^{7}$ | $4.36\cdot 10^{7}$
$24\times 14$ | $1978$ | $182$ | $43236$ | $1.05$ | $8.83\cdot 10^{7}$ | $1.05\cdot 10^{8}$

The algorithmic complexity heavily depends on the implementation, computer
architecture, and type of fitness function. In particular, the minimization of
Q-factor requires the evaluation of quadratic forms representing stored energy
and reactive power, which is computationally expensive. From this perspective,
the minimization of Q-factor represents a computationally challenging example.

Figure 3: Performance of the global step for Q-factor minimization from Fig.
1. The mesh grid of $12\times 7$ pixels (dashed lines in Fig. 1) is used for
electrical size $ka=0.5$. The performance is normalized to the fundamental
bound $Q_{\mathrm{lb}}^{\mathrm{TM}}$. Three implementations are compared:
memetic combination of global and local steps proposed in this paper (red
curves), the sole global step (a sole heuristic genetic algorithm, blue
curves) and adjoint formulation of density-based topology optimization (black
curves). Each marker represents one global iteration of the optimizer. All
algorithms were run five times. The filled areas (yellow, orange, and violet)
show the period when the local step of memetic scheme is performed in
$j=\left\\{2,3,5\right\\}$ global iteration, cf., Fig. 4.

A graphical representation of the performance of the memetic algorithm (global
and local steps together) is shown in Fig. 3 for a $12\times 7$ mesh grid and
electrical size $ka=0.5$, with computational time used as the $x$-axis. Only
$N_{\mathrm{A}}=16$ agents are required for the memetic algorithm to reach
very low values of the optimized metric. The remarkable convergence rate of
the memetics stems from the local step. The local step, when it is applied for
the first time in the second global iteration ($j=2$), already finds a
solution (after $30$ s) close to the final result. This occurs because of
properties of the Q-factor’s solution space, which makes it possible to find a
reasonable solution from an arbitrary starting point, see Monte Carlo analysis
in [14].

Figure 4: Performance of the local step for three different global step
iterations $j=\left\\{2,3,5\right\\}$. A mesh grid with $12\times 7$ pixels,
i.e., the same as in Fig. 3, is used. The performance is normalized to the
fundamental bound $Q_{\mathrm{lb}}^{\mathrm{TM}}$. All $N_{\mathrm{A}}=16$
agents are shown in each global iteration $j$. The computational time spent
with the local updates in global iteration $j$ is highlighted by the same
background color (yellow, orange, violet) in this figure and in Fig. 3. The
thick red curve represents an average of all $N_{\mathrm{A}}$ agents (only for
computational times $t$ where all curves exist). Each marker represents one
local iteration of the optimizer. The minimum, maximum, and average values are
shown for the averaged cost function, their beginning and end.

The power of the local step within the memetic scheme is emphasized in Fig. 4
where the cost function for all $N_{\mathrm{A}}=16$ agents of memetics from
Fig. 3 are depicted for global iterations $j\in\left\\{2,3,5\right\\}$.
Iteration $j=2$ (the top pane) is the first one when the local step is used
(the first iteration only evaluates random initial seeds, see Part I [11, Fig.
8]). Although it takes the majority of the computational time, it is also
capable of decreasing the value of the objective function by three orders in
magnitude, reaching, in one case, $Q/Q_{\mathrm{lb}}^{\mathrm{TM}}\approx
1.15$. Then the next iteration (the middle pane) starts, on average, at a
higher value of the objective function since the global operators were applied
over the resulting words, consequently perturbing them and increasing the
chance that a better local optimal will be found. Subsequently, the objective
function is improved by one to two orders at the expense of a few seconds.
Later, the fifth iteration slightly improves the local minima, almost reaching
the final values of the objective function found for this grid, cf., Table I.

###  III-D Comparison with Other Optimization Schemes

Apart from the performance of the memetic algorithm, the results obtained by
sole heuristics (pure genetic algorithm, only the global step) and by density-
based topology optimization [2] are also presented in Fig. 3. An amount of
$N_{\mathrm{A}}=80$ agents were required for the genetic algorithm to reach
values comparable to those offered by the memetic algorithm. The density-based
topology optimization was using the method of moving asymptotes [42] with the
interpolation function from [7]. The density filter [43] had a fixed radius of
a linearly decaying cone. A continuation scheme [44] was utilized to ensure
convergence to the structure consisting of either PEC or vacuum only.

The performance offered by the sheer genetic algorithm (blue lines) is
inferior to the memetic scheme.

The classical topology optimization, represented by black lines in Fig. 3,
excels in speed. Consequently, it will supersede the memetic scheme when the
number of DOFs is scaled up. This can be mitigated by using techniques such as
surrogate modeling [45, Fig. 2]. Density-based topology optimization is,
however, not able to reach as low values of the optimized metric as compared
to the memetic scheme. This is mostly given by the necessity of thresholding
that converts the resulting “gray” design to a binary design made solely of
PEC and vacuum. The abrupt jumps in Fig. 3 are caused by changing the
properties of the filter in order to gradually remove the “gray”-scaled
elements at the course of the optimization [2]. Another drawback of classical
topology optimization is the necessity to know the proper shape of the density
curve and the size and shape of density filters. These are heuristic
parameters not present in the memetic scheme. A notable strength of the
memetic algorithm is that explicit knowledge of derivatives of the optimized
metric with respect to material variables is not needed.

###  III-E Spherical Shell as a Design Region

Spherical helices fed by a delta gap source can simultaneously excite both TM
and TE modes [46, 37], thereby reaching the lower bound $Q_{\mathrm{lb}}$ [47,
39] and breaking the $Q_{\mathrm{lb}}^{\mathrm{TM}}$ bound. To construct them
ad hoc without a long study of the topic and deep empirical knowledge is,
nevertheless, difficult [48]. In this section, the Q-factor without the self-
resonance constraint is minimized for a spherical shell of electrical size
$ka=0.2$, and both $Q_{\mathrm{lb}}$ and $Q_{\mathrm{lb}}^{\mathrm{TM}}$
bounds222Considering electrically small region, $ka\ll 1$, the bounds for a
spherical shell of radius $a$ are $(ka)^{3}Q_{\mathrm{lb}}\approx 1$ and
$(ka)^{3}Q_{\mathrm{lb}}^{\mathrm{TM}}\approx 3/2$ [49]. are utilized to judge
the performance of the optimized structures.

The spherical shell was discretized into $N=2304$ DOF and optimized with
$N_{\mathrm{A}}=64$ agents and $J=250$ global iterations. The relative
differences, both for global and local steps, were set once again to
$\epsilon=10^{-7}$. One discrete delta gap source was used (its position and
orientation is irrelevant thanks to the spherical symmetry); the materials are
either PEC or vacuum.

The optimization ran for $46$ hours and the resulting structure is shown in
Fig. 5. While the computational time might seem enormous, in total $6.24\cdot
10^{8}$ antennas were evaluated via the exact reanalysis local step (within a
solution space containing $2^{2303}\approx 1.87\cdot 10^{693}$ antenna
variants). The crucial role of the local step is, once again, visible in Fig.
6 where the first iteration with a greedy search decreases the objective
function from $(ka)^{3}Q\approx 10^{3}$ to $(ka)^{3}Q\approx 1.4$. Good
performance in the Q-factor, similar to what was reported in [37], where the
antenna was designed empirically, is observed, i.e., $Q/Q_{\mathrm{lb}}\approx
1.27$ and $Q/Q_{\mathrm{lb}}\approx 0.91$. The helix has a slope optimized for
an optimal combination of TM and TE modes [39] and several turns to reach
self-resonance. The input impedance is, however, only $Z_{\mathrm{in}}\approx
2.14\,\Omega$, which seriously limits the practical usage of such an antenna.
It can be increased by prioritizing low reflectance instead of Q-factor in a
multi-objective optimization as shown in the next section.

Figure 5: (left) A structure optimized within a design region coinciding with
a spherical shell of electrical size $ka=0.2$, $N=2304$. The structure is fed
by a delta gap feeder depicted by the thick blue line and it resembles a
spherical helix, cf. [46]. Normalized Q-factor reaches $(ka)^{3}Q\approx
1.05$, i.e., a value below the TM-bound
($Q/Q_{\mathrm{lb}}^{\mathrm{TM}}\approx 0.91$) and slightly above the
fundamental bound ($Q/Q_{\mathrm{lb}}\approx 1.27$). (right) Current density
on the structure from the left pane. The absolute value of the current density
is depicted by the colormap, the direction of the real part of the current is
depicted by arrows. These plotting settings are used throughout the paper.
Figure 6: The cost function of the first ten selected agents (out of $64$)
based on net computational time (in hours). Due to major changes in the first
iterations, two panes of different $y$-axis limits are shown. The cost
functions do not decrease monotonically during the global step since
heuristics perturb the reached local minima. The thick red curve shows the
minimum reached in each iteration of global step $j$. The inset in the top-
left corner shows the design region fully populated with DOF at the beginning
of the optimization, and the bottom right inset shows the optimized structure,
cf., 5.

##  IV Trade-off Between Q-factor and Input Impedance

It was demonstrated in the previous example that more than one metric is of
concern in antenna design. For example, to minimize the Q-factor and to
reasonably match the antenna to a given input impedance
$Z_{\mathrm{in}}=R_{\mathrm{in}}+\mathrm{j}X_{\mathrm{in}}$ simultaneously
requires composite objective functions. It can, for example, be of the
following form

$f\left(\mathbf{g}\right)=\dfrac{Q\left(\mathbf{g}\right)}{Q_{\mathrm{lb}}}\left(1+\zeta\left|\varGamma\left(\mathbf{g},Z_{\mathrm{in}}^{0}\right)\right|^{2}\right),$
(3)

where $\varGamma\left(\mathbf{g},Z_{\mathrm{in}}^{0}\right)$ is the reflection
coefficient

$\varGamma\left(\mathbf{g},Z_{\mathrm{in}}^{0}\right)=\dfrac{Z_{\mathrm{in}}(\mathbf{g})-Z_{\mathrm{in}}^{0}}{Z_{\mathrm{in}}(\mathbf{g})+Z_{\mathrm{in}}^{0}}$
(4)

of a structure represented by a word $\mathbf{g}$, evaluated for
characteristic impedance $Z_{\mathrm{in}}^{0}$ [50], and $\zeta$ is a
weighting coefficient. For $\zeta=0$ the optimization is the same as in
Section III. Any value $\zeta>0$ takes into account the matching as well. This
results in a Pareto-type optimization with a scalarized objective function
[51], where the parameter $\zeta$ sweeps over the feasible set.

The PEC rectangular plate of electrical size $ka=0.7$ from Section III-B is
considered first. Self-resonance is reachable for this electrical size, and
only the real part of the input impedance has to be optimized. The desired
characteristic impedance is set to $Z_{\mathrm{in}}^{0}=50\,\Omega$ and weight
$\zeta$ is swept from $\zeta=0$ to $\zeta=5$ in $41$ equidistant samples. The
results are depicted in Fig. 7. The Pareto-optimal solutions [52] are
highlighted by red circles and interconnected to form an approximate Pareto
frontier. It is seen that the Q-factor, in terms of
$Q/Q_{\mathrm{lb}}^{\mathrm{TM}}$, can be minimized to
$Q/Q_{\mathrm{lb}}^{\mathrm{TM}}\approx 1.02$ but at that point the antenna is
not matched. Conversely, a slight increase in
$Q/Q_{\mathrm{lb}}^{\mathrm{TM}}$ leads to excellent matching. The two most
distinct solutions, marked by (A) and (B) in Fig. 7 are shown in Fig. 8 in
terms of optimal structure (top) and current density (bottom). The left
structure in Fig. 7, denoted as (A), performs best in terms of Q-factor. The
right structure performs well in Q-factor and is matched to $50\,\Omega$. The
visual comparison reveals that the parallel stub was created in structure (B)
to match the input impedance. This technique [36] is also used in practice
[37, Fig. 9].

Figure 7: Trade-off between low Q-factor $Q/Q_{\mathrm{lb}}^{\mathrm{TM}}$ and
matching to $Z_{\mathrm{in}}^{0}=50\,\Omega$. A PEC rectangular plate of
$12\times 7$ pixels, discretized into $N=485$ DOF was used. The electrical
size is $ka=0.7$. The delta gap feeder is placed in the top middle. In total,
$41$ samples were evaluated with weight $\zeta$ in (3) set from 0 to 5 with an
equidistant step. The Pareto frontier is highlighted by the red dashed line.
Pareto-optimal solutions are shown as red marks. The two extreme cases, the
one with the lowest Q-factor and the one well-matched to $Z_{\mathrm{in}}^{0}$
are shown in Fig. 8. Figure 8: The optimal shapes (top) and surface current
densities (bottom) for two members of the Pareto frontier in Fig. 7. Enabled
DOF and material used are highlighted. The structures are fed by a delta gap
feeder, depicted by the thick blue line. The red lines denote enabled DOF. The
(A) case minimizes Q-factor and resembles a meanderline. The (B) case is
matched, therefore, there is a stub close to the feeder.

The same formulation (3) is applied to an optimization of the spherical shell
presented in Section III-E with an attempt to modify the structures from Fig.
5 to be matched to $Z_{\mathrm{in}}^{0}=50\,\Omega$. For this purpose, weight
$\zeta$ in (3) is set to $\zeta=1$. The optimization found a minimum with
$Q/Q_{\mathrm{lb}}\approx 1.27$ and $Q/Q_{\mathrm{lb}}^{\mathrm{TM}}\approx
0.91$, i.e., the same as before, and
$Z_{\mathrm{in}}\approx(51.5-\mathrm{j}2.39)\,\Omega$, i.e.,
$|\varGamma|^{2}\approx 7.72\cdot 10^{-4}$, which is sufficient for a majority
of applications [53]. The optimal structure is shown in Fig. 9. When the
optimal shapes and currents for an unmatched antenna in Fig. 5 and a matched
antenna in Fig. 9 are compared, a similar difference, as in Fig. 8, is
observed, i.e., the optimal structure is only slightly perturbed in the
vicinity of the feeder, using a stub-matching technique to divide the current
flowing through the feeding port. Another improvement is a wider metallic
strip used to curve the helix antenna (at least two DOF are enabled per width
of the strip).

Figure 9: (left) A structure optimized within a design region coinciding with
a spherical shell of electrical size $ka=0.2$, $N=2304$. The structure is fed
by a delta gap depicted by the thick blue line. The optimization considers
Q-factor minimization and matching to $Z_{\mathrm{in}}^{0}=50\,\Omega$.
(right) Current density on the structure from the left pane. Similar to Fig.
5, normalized Q-factor reaches $(ka)^{3}Q\approx 1.05$, i.e.,
$Q/Q_{\mathrm{lb}}\approx 1.27$, and $Q/Q_{\mathrm{lb}}^{\mathrm{TM}}\approx
0.91$, but the reached input impedance is
$Z_{\mathrm{in}}\approx(51.5-2.39\mathrm{j})\,\Omega$ ($|\varGamma|^{2}\approx
7.83\cdot 10^{-4}$).

Considering Q-factor, this section shows that there are many local minima
close to the global minimum. However, other antenna metrics, such as input
impedance, have different values in these local minima. In other words,
Q-factor minimization can always be extended to a multi-objective case with a
constraint on another metric without suffering a major increase in Q-factor.
This is consistent with empirical evidence known in the literature [37].

##  V Maximal Realized Gain of an Array

Electrically larger structures are considered in this example. For this
reason, the structure is preselected (an antenna array) and optimized further.
A thin-strip linear antenna array consisting of $N_{\mathrm{dip}}$ dipoles,
operating at $1\,$GHz, and made of copper, $\sigma=5.96\cdot
10^{7}\,\mathrm{Sm}^{-1}$, is generated. The thin-sheet model is utilized333It
is assumed that the thickness of the structure is much larger than the
penetration depth. The current is modeled as exponentially decaying from the
surface of infinite half-space [54]. to determine the surface resistivity as

$R_{\mathrm{s}}=\sqrt{\dfrac{kZ_{0}}{2\sigma}},$ (5)

where $Z_{0}$ is the impedance of free space. The initial length of the
dipoles is $\ell/\lambda=0.55$ and is intentionally longer than
$\ell=\lambda/2$ to give the algorithm a chance to adopt it freely. The width
of the strips is $\ell/60$, and the separation distance is $d$. The first
element from the left always acts as a passive reflector, while the second
dipole from the left is always fed by the delta gap feeder in the middle.
Other dipoles act as passive directors.

Maximal realized gain $G_{\mathrm{r}}$ is chosen as a figure of merit. Since
the optimization algorithm performs the minimization, the sign is switched to
minus as

$f\left(\mathbf{g}\right)=-G_{\mathrm{r}}\left(\mathbf{g}\right)=-G(\mathbf{g},\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})\left(1-\left|\varGamma\left(\mathbf{g},Z_{\mathrm{in}}^{0}\right)\right|^{2}\right)$
(6)

where $G(\mathbf{g},\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})$ is antenna
gain

$G(\mathbf{g},\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})=\dfrac{U(\mathbf{g},\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})}{P_{\mathrm{rad}}(\mathbf{g})+P_{\mathrm{lost}}(\mathbf{g})}=\dfrac{U(\mathbf{g},\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})}{P_{\mathrm{tot}}(\mathbf{g})}$
(7)

and $\varGamma\left(\mathbf{g},Z_{\mathrm{in}}^{0}\right)$ is reflection
coefficient (4). From an implementation point of view, both radiation
intensity and total power are one order cheaper to evaluate than Q-factor. As
the strip (1D) structures are optimized, we can, therefore, expect relatively
fast convergence.

Figure 10: Shape optimization of an antenna array consisting of
$N_{\mathrm{dip}}=\left\\{3,4,5\right\\}$ thin-strip dipoles, made of copper
($\sigma=5.96\cdot 10^{7}\,\mathrm{Sm}^{-1}$) and operating at $1$ GHz. The
initial length of the optimized dipoles was $\ell/\lambda=0.55$ and the width
$w=\ell/60$. The optimization was repeated for various separation distances
denoted as $d/\lambda$ to maximize realized gain. The reference impedance was
set to $Z_{\mathrm{in}}^{0}=50\,\Omega$, polarization pointed along the
dipoles, and gain was measured in the end-fire direction. Excitation was
performed by a delta gap feeder placed in the middle of the second dipole from
the left, see Fig. 11. The vertical dotted lines point to the maximal realized
gain found for a various number of dipoles and denote the corresponding
electrical length $kd$. For convenience, the results are compared with
fundamental bounds with prescribed input impedance (dashed lines). For one
optimized sample, see Fig. 11.

The minimization of (6) was repeated for
$N_{\mathrm{dip}}=\left\\{3,4,5\right\\}$ dipoles with separation distance
selected from $d/\lambda=0.02$ to $d/\lambda=0.5$, see the solid lines in Fig.
10. To be as thorough as possible, the fundamental bound on the realized gain
with $Z_{\mathrm{in}}^{0}=50\,\Omega$ was also evaluated (see the dashed
lines). The procedure is detailed in Appendix B.

There are many conclusions that can be drawn from Fig. 10. First, there is a
“sweet spot” in separation distance $d/\lambda$ (or in $kd$ which denotes the
electrical length of an array) and this optimal separation distance is
different for optimized structures and for fundamental bounds (fundamental
bounds reach maxima for lower distance $kd$). Second, a higher number of array
elements leads to significantly higher realized gain, see Table II for
numerical comparison. This is consistent with array theory [55]. Finally, the
maxima given by the fundamental bounds are closely followed by realized
designs. As the desired input impedance was prescribed for the bounds, see
Appendix B, this might indicate that all the realized arrays were sufficiently
well matched and that the antenna gain for both bounds and designs are
comparable. This has been verified by an inspection of the optimization data.

TABLE II: Comparison of realized gain maxima found via memetic shape optimization and the fundamental bounds depending on the number of dipoles used. Both separation distance $d/\lambda$ and electrical length of the array $kd$ are shown | memetics | fundamental bound
---|---|---
$N_{\mathrm{dip}}$ | $d/\lambda$ | $kd$ | $G_{\mathrm{r}}$ | $d/\lambda$ | $kd$ | $G_{\mathrm{r}}$
$3$ | $0.260$ | $3.27$ | $\mathbf{8.13}$ | $0.200$ | $2.51$ | $\mathbf{10.4}$
$4$ | $0.290$ | $5.47$ | $\mathbf{12.0}$ | $0.245$ | $4.62$ | $\mathbf{15.6}$
$5$ | $0.320$ | $8.04$ | $\mathbf{15.7}$ | $0.275$ | $6.91$ | $\mathbf{21.1}$

As an example, the optimal design for $N_{\mathrm{dip}}=4$ dipoles and a
separation distance chosen for the highest realized gain ($d/\lambda=0.29$) is
shown in Fig. 11. It is seen that the structure was significantly modified. A
driven element was split by the removal of one DOF. This is indicated by the
gray dashed line and a tag “cut” in Fig. 11. The same result occurred for the
first and fourth elements (counted from the left). The third element was
modified by the complete removal of material from its lower part. It is
obvious that these modifications cannot be achieved by a simple parametric
sweep to optimize the overall length of the dipoles. Such an approach would
remove all material below the “cut” label, reducing the performance from
$G_{\mathrm{r}}\approx 12.0$ to $G_{\mathrm{r}}\approx 10.9$. The resulting
structure is not symmetrical, even though the initial problem was. This is
partly because the underlying discretization grid is not symmetrical and
partly because the optimization problem is non-convex.

Figure 11: Optimized array of four thin-strip dipole elements made of copper
for maximal realized gain at $1$ GHz. The separation distance was fixed to
$d/\lambda=0.29$, i.e., for a distance where the highest realized gain is
reached, cf., Fig. 10. The array is fed in the middle of the second dipole to
the left. The reference impedance is $Z_{\mathrm{in}}^{0}=50\,\Omega$,
polarization $\hat{\boldsymbol{e}}=\hat{\boldsymbol{z}}$, and direction of
main lobe $\hat{\boldsymbol{d}}=\hat{\boldsymbol{x}}$. The initial electric
length $\ell/\lambda=0.55$ is modified by cutting three dipoles. One dipole is
modified by the complete removal of the material.

##  VI Maximum Absorption in a Given Region

The last example deals with the maximization of power absorbed in a given
region, e.g., an optimal receiving structure for an RFID chip. The
optimization domain is, in this case, only a part of the entire structure, see
Fig. 12 for a schematic layout. The controllable part, which is a subject of
the optimization (highlighted by the yellow color) is made of copper,
$\sigma=5.96\cdot 10^{7}\,\mathrm{Sm}^{-1}$. The chip (pink color) is made of
carbon, $\sigma=1\cdot 10^{4}\,\mathrm{Sm}^{-1}$. The thin-sheet resistivity
model (5) is applied. The operational frequency is $3\,$GHz, the shorter side
has length $\lambda/4$ and the longer $11\lambda/32$, and the incident plane
wave is impinging from the end-fire direction, i.e., along the
$-\hat{\boldsymbol{x}}$ direction with
$\hat{\boldsymbol{e}}=\hat{\boldsymbol{y}}$ polarization, see Fig. 12.

Figure 12: A schematic for maximization of power absorbed in the prescribed
region made of carbon ($\sigma=1\cdot 10^{4}\,\mathrm{Sm}^{-1}$). The
operation frequency is $f=3\,$GHz, the shorter side has length $\lambda/4$.
The amount of power absorbed from the plane wave impinging from the
$\hat{\boldsymbol{d}}=-\hat{\boldsymbol{x}}$ direction with
$\hat{\boldsymbol{e}}=\hat{\boldsymbol{y}}$ polarization is evaluated only in
this region denoted by the pink color. The optimization domain is made of
copper ($\sigma=5.96\cdot 10^{7}\,\mathrm{Sm}^{-1}$) and denoted by the yellow
color.

The absorbed power $P_{\mathrm{lost}}$ is evaluated from

$P_{\mathrm{lost}}=\dfrac{1}{2}\mathbf{I}^{\mathrm{H}}\left(\mathbf{D}_{\mathrm{chip}}^{\mathrm{H}}\mathbf{R}_{\rho}\mathbf{D}_{\mathrm{chip}}\right)\mathbf{I},$
(8)

where $\mathbf{I}$ is the current on the optimized structure,
$\mathbf{R}_{\rho}$ is the lossy matrix, and matrix
$\mathbf{D}_{\mathrm{chip}}$ is an indexation matrix having zeros everywhere
except for diagonal positions corresponding to the DOF lying in the fixed
region, see Fig. 12. To ease the computational burden, the matrix
$\mathbf{D}_{\mathrm{chip}}$ is used to index out only the relevant entries
required to evaluate (8).

The upper bound for this optimization problem is
$P_{\mathrm{lost}}^{\mathrm{ub}}\approx 11.2$ µW and was found using a
procedure specified in Appendix C. The optimal current impressed in vacuum is
depicted in Fig. 13. Even though this current does not fulfill
$\mathbf{Z}\mathbf{I}=\mathbf{V}$ with a plane wave excitation represented by
a vector of expansion coefficients $\mathbf{V}$, it is seen that it is a
smooth function with maxima along the shorter sides and in the area of the
chip. It must be noted that there is no guarantee that the bound is tight in
this case.

Figure 13: Current density representing the optimal performance for maximal
power absorbed $P_{\mathrm{lost}}$ in the uncontrollable region, depicted by
the pink color in Fig. 12.

The optimization was performed twice using the same settings and only varying
the number of agents $N_{\mathrm{A}}$ used for the global step. As in all
previous examples, the first application of the local step led to an immense
improvement of the objective function value, here, of absorbed power
$P_{\mathrm{lost}}$ from approximately $3.50\cdot 10^{-3}$ µW to $5.10$ µW,
i.e., by more than three orders in magnitude. For $N_{\mathrm{A}}=80$ agents
(five times the number of cores of the Threadripper 1950X), the optimization
ran for approximately $4.3$ hours with maximum $P_{\mathrm{lost}}\approx 6.96$
µW, see Fig. 14. Increasing the number of agents to $N_{\mathrm{A}}=192$ led
to maximum $P_{\mathrm{lost}}\approx 7.75$ µW found in $10.9$ hours. This
value reaches $67.6\%$ of the upper bound realized by the current shown in
Fig. 13. As compared to the original structure depicted in Fig. 12, the power
absorbed in the “chip” was increased by a factor of $4.52\cdot 10^{4}$.

Figure 14: Cost function for maximization of absorbed power
$P_{\mathrm{lost}}$. Optimization was performed for $N_{\mathrm{A}}=80$ (blue
line) and $N_{\mathrm{A}}=192$ (red line) agents. The optimal structures for
both cases are shown in the bottom right corner. The structure for
$N_{\mathrm{A}}=192$ is shown in detail in Fig. 15.

Comparing the objective functions in Fig. 14, it is clear that increasing
$N_{\mathrm{A}}$ increases the number of simultaneously used local minima
during the evaluation of the global step, which leads to increased diversity
and, consequently, increases the chance of finding a high-quality solution.
This is, of course, at the expense of computational time, where the increase
is approximately linear with $N_{\mathrm{A}}$. Notice, however, that with
access to the computational cluster with $C$ cores, there is almost no
difference between computational times as far as $N_{\mathrm{A}}\leq C$
because all agents in the same global iteration $j$ are evaluated
simultaneously thanks to excellent scalability in parallel computing.

The optimal candidates for both runs are depicted in the bottom-right corner
of Fig. 14 which shows that the shape for $N_{\mathrm{A}}=192$ is more
regular, consisting of three wire structures. The shape for
$N_{\mathrm{A}}=80$ has not only poorer performance in terms of absorbed power
$P_{\mathrm{lost}}$, but is also significantly more complex rendering its
possible manufacturing risky.

Figure 15: (left) A structure optimized within a design region depicted in
yellow in Fig. 12 for the maximization of power absorbed in a chip depicted in
pink in Fig. 12. The physical setting is as described in Fig. 12. The number
of DOF is $N=905$, the number of controllable DOF is $B=863$. The number of
agents used was $N_{\mathrm{A}}=192$, see the red curve in Fig. 14. The red
connections denote enabled DOF. (right) Current density on the optimized
structure from the left pane.

##  VII Discussion and Future Challenges

The method has some unique properties and offers unorthodox features which are
discussed below together with a list of pros and cons.

###  VII-A Properties and Features

####  VII-A1 Full-wave evaluation

The entire approach is full-wave. The numerical errors occurring during the
iterative updates are negligible, typically influencing the last digit in
double precision444This was verified by optimizing the shape with the proposed
method first and solving the method of moments for the resulting shape
directly after that. The difference in objective function was compared then
for these two current densities..

####  VII-A2 Fixed discretization grid

The procedure is based on a fixed discretization grid. It may happen that some
objective function tends to thin structures, i.e., only one DOF per width is
chosen no matter what density of grid is used. A procedure with remeshing can
be, however, applied in these cases. The objective function can also be
penalized with ohmic losses, the incorporation of which always reflects
physical reality better.

####  VII-A3 No interpolation function

As compared to the classic topology optimization [2, 8], there is no
interpolation procedure which can change the actual value of an objective
function when performed and which selection depends on a user.

####  VII-A4 Gradient-based procedure

As opposed to the pixeling [56], the proposed method involves the local step.
Gradient-based optimization methods are typically preferred due to their rapid
convergence [57]. The gradient-based method may fail if a problem includes a
significant amount of local minima with the solution space unable to be
regularized [20]. To mitigate this, the global step is utilized. Being
connected is, however, a question of whether the greedy search used in the
local step is an appropriate choice as there are NP-hard problems known to be
“greedy-resistant” [58]. Here we only rely on the experimental evidence
presented in this paper, which suggests that the greedy search provides
appropriate designs.

####  VII-A5 Second topology derivative

Considering genetic algorithm as a particular choice for the global step, the
mutation operator [59] serves as a second option. It is performed for random
perturbations, but its application can be controlled by a proper choice of
metaparameters.

####  VII-A6 MoM compatibility

The approach is compatible with any MoM formulation based on piecewise basis
functions and a solution found with a direct (noniterative) solver. Important
differences may, however, exist depending on the type of basis functions used.
For example, overlapping basis functions (e.g., Rao-Wilton-Glisson (RWG)
[26]), do not coincide one by one with discretization elements, therefore,
removing one DOF does not correspond to removing one discretization element.

####  VII-A7 Compatibility with Fundamental Bounds

The formulation within the framework of MoM allows evaluating fundamental
bound on many metrics of practical interest. These bounds then serve as
stopping criteria for topology optimization, a non-convex problem in which the
proximity to the global extremum is unknown.

####  VII-A8 Multi-modality

The algorithm is multi-modal since, in principle, $N_{\mathrm{A}}$ local
minima are found during each global step thanks to the application of the
local step [60]. The number of unique local minima decreases as the
optimization converges to the global minimum.

####  VII-A9 Big data

A huge amount of data is gathered during the optimization, cf. Table I which,
together with excellent control over the entire optimization, makes this
approach an ideal candidate for real-time data processing with machine
learning.

####  VII-A10 Flexible objective function

Topology sensitivity is evaluated as the difference between the performance of
the actual current $\mathbf{I}\left(\mathbf{g}_{i}\right)$ and the performance
of the current flowing on the perturbed structure
$\mathbf{I}\left(\mathbf{g}_{i+1}\right)$. As such, the objective function can
be easily extended towards:

  * •

an objective function evaluated only within a sub-region (see Section VI),

  * •

an objective function taking into account multiple frequency samples,

  * •

an objective function involving port quantities defined for multiple ports.

####  VII-A11 Multi-objective optimization

Multi-objective optimization can easily be performed with the proposed
memetics by scalarization, see Section. IV. Only convex trade-offs can,
however, be found with a weighted sum, such as in (1) or (3) where more
advanced techniques would have to be adopted for complex Pareto frontiers,
e.g., a rotated weighted metric method [60].

###  VII-B Generalization of Variable Space

####  VII-B1 Feeding synthesis

The optimal placement of the feeder(s), including optimal amplitudes and
phases, can be incorporated into shape optimization. For one feeder, it is,
technically, not needed555It is assumed that the optimal shape is adjusted
according to the initial placement of the feeder.. For more feeders, it is an
extra combinatorial problem which can run over or run together with the shape
optimization.

####  VII-B2 Optimization variables

It is possible to generalize the entire algorithm so that discretization
elements are dealt with and not DOF (basis functions). One possibility is to
define the smallest blocks manually, and optimize over them only. This slows
down the evaluation since vectorization is not as efficient because block
inversion has to be applied. On the other hand, the number of unknowns is
reduced depending on the number and size of the optimized blocks. The initial
tests indicate that computational time is comparable or lower than with DOF as
the optimization unknowns.

####  VII-B3 Optimization domain

Topology sensitivity can be evaluated only for removals or additions offering
a substantial speed-up. Alternatively, as often required in practice, only a
part of the structure can be optimized, see Fig. 12, reducing computational
time proportionally to the size of the optimized region. Another possibility
is to investigate only DOF, which separates two regions filled by different
materials (the number of unknowns is reduced approximately by one order).

####  VII-B4 Multi-material representation

The two-state, binary, material representation (PEC $\times$ vacuum) can
easily be generalized to a multi-state optimization at the cost of the linear
increase of variable space size.

###  VII-C Deficiencies and Possible Remedies

####  VII-C1 Slow evolution

Shape modification is performed via the removal or addition of a DOF
decreasing the value of the objective function the most. Following this greedy
search, the DOF corresponding to the thickest green line in Fig. 16, i.e., the
connections located to the left and on the top right section of the feeder are
to be removed in order to eliminate the short-circuit. Such an approach is,
however, relatively slow666As compared to the adjoint formulation of topology
optimization [2] where all DOF are updated at once., requiring in many cases
${\mathcal{O}}(N)$ local updates. The possible remedy might be to average the
sensitivity in a small area and update all DOF inside at once.

Figure 16: Topology sensitivity map of a rectangular PEC plate of electrical
size $ka=0.5$ and a hole in the middle. The objective function evaluates
Q-factor (1). The delta gap feeder is depicted by the thick blue line. The
sensitivites $\tau_{n}$ are evaluated for all $N$ DOF. The negative
sensitivities, $\tau_{n}<0$, are highlighted by the green lines, the positive
sensitivies, $\tau_{n}>0$, by the red lines. The thicker the line, the higher
the absolute value of the sensitivity.

####  VII-C2 Algorithm complexity

Numerical complexity is unpleasant as it favors small to mid-size structures.
Realistically, with the current implementation and state-of-the-art hardware,
thousands of unknowns are possibly optimized in hours. Billions of DOF are,
however, common in FEM [9], but this is because of different algebraic
properties of the stiffness matrix as compared to the impedance matrix.

####  VII-C3 Shape irregularity

The method might produce highly irregular optimized shapes. This is a common
issue of shape optimization, no matter what technique is used, see Fig. 17,
comparing representatives of this method, pixeling based on genetic algorithm
[56], and gradient-based topology optimization [8]. All shapes are relatively
irregular, i.e., this property is common to shape optimization routines.
Hence, the next step could be to find a way how to regularize the shapes.

Figure 17: Comparison of final structures optimized with various shape
optimization techniques. (a) The method proposed in this work, Q-factor
minimization. (b) Pixeling based on the genetic algorithm, Q-factor
minimization [61]. (c) Adjoint formulation of topology optimization, total
efficiency minimization [8]. All techniques utilize method of moments to
evaluate the objective function. The orange dashed line delimits the design
region.

##  VIII Conclusion

A novel memetic procedure for optimizing electromagnetic devices was
presented, consisting of local and global steps to mitigate the disadvantages
of both approaches when used alone. The fixed discretization grid is assumed,
allowing the method of moments system matrix to be inverted only once while
storing it and other required matrices in the computer’s memory. The iterative
full-wave evaluation of all the smallest topology perturbations is done using
an inversion-free (Sherman-Morrison-Woodbury) formula and an exact-reanalysis-
based procedure. The local step provides gradient-type information about the
topology, and local updates are performed in a greedy sense, i.e. with the
part of the optimized shape enhancing the performance the most being updated
in each local iteration. The global step used in this paper is based on a
genetic algorithm which maintains diversity and increases the chance of the
algorithm’s convergence towards a global minimum. Genetic algorithms match
well with the discrete form of the optimization problem and, as a result, the
memetic procedure is robust, fast, and versatile.

The procedure has many advantages and unique features, such as being able to
find minima close to the fundamental bounds, implying that the local minima
are close to the global one. The number of full-wave evaluated shapes spans
from millions for small problems to billions for larger problems with
thousands of unknowns. The optimization approach is multi-modal and capable of
finding many local minima at once during one global step. As a result, the
global scheme optimizes shapes only in a significantly reduced solution space
containing only local minima. Since the optimized shape is known at every
step, a multi-port, multi-frequency, or multi-material optimization can be
performed only at the expense of the linear increase of computational time.

Examples shown in this part proved the efficiency of the method. For the first
time, the performance of optimized structures has been directly compared with
the fundamental bounds on the optimized metrics. This practice provides an
ultimate measure of optimization efficiency, but also naturally scales
optimized metrics and establishes the straightforward terminal criterion for
optimization. In many cases, the performance of the resulting shapes closely
follows the bounds. The excellent performance of the algorithm demonstrates
its wide applicability in various inverse design problems in electromagnetics.

There are still many future challenges that could broaden the usage of the
proposed method to make it even more efficient, such as including lumped
elements and their synthesis. It can be shown that this technique is similar
to adjoint formulation over gray-scaled material. Another possibility is to
introduce a set of geometrical operators, representing geometrical metrics,
such as the area spanned by the material, or the curvature of the shape.
Having complete freedom in formulating the objective function, these operators
may improve the regularity of shapes and eliminate manufacturing difficulties.
Technical, and still extremely relevant, improvements include GPU
implementation or a detailed study of metaparameters used for optimization
settings and their tuning at the beginning or throughout an optimization.
Another way to accelerate the evaluation is to use an adaptive scheme with a
successively refined discretization grid to impose the optimal structure from
a coarse grid as the initial shape for the finer grid.

##  Appendix A Fundamental Bound on Antenna Q-factor

The lower bound on radiation Q-factor $Q_{\mathrm{lb}}$ is found via a
quadratically constrained quadratic program (QCQP)

$\begin{split}\underset{\mathbf{I}}{\mathrm{minimize}}\quad&\mathbf{I}^{\mathrm{H}}\mathbf{W}\mathbf{I}\\\
\mathrm{subject\,to}\quad&\mathbf{I}^{\mathrm{H}}\mathbf{R}_{0}\mathbf{I}=1\\\
&\mathbf{I}^{\mathrm{H}}\mathbf{X}\mathbf{I}=0,\end{split}$ (9)

where $\mathbf{R}_{0}$ is radiation and $\mathbf{X}$ is the reactance part of
impedance matrix $\mathbf{Z}$, and $\mathbf{W}$ is the stored energy matrix.
The problem is recast into its dual form and solved via a generalized
eigenvalue problem as described in [62].

Q-factor $Q_{\mathrm{lb}}^{\mathrm{TM}}$ is a tighter bound for all cases when
only the TM modes are involved. This includes, e.g., planar antennas with a
discrete feeder [35, 31], i.e., the case studied in Section III. The formula
(9) still applies with a change
$\mathbf{R}_{0}\to\mathbf{R}_{0}^{\mathrm{TM}}$, where matrix
$\mathbf{R}_{0}^{\mathrm{TM}}$ is evaluated as
$\mathbf{R}_{0}^{\mathrm{TM}}=\left(\mathbf{U}_{1}^{\mathrm{TM}}\right)^{\mathrm{H}}\mathbf{U}_{1}^{\mathrm{TM}}$
with $\mathbf{U}_{1}^{\mathrm{TM}}$ being a projection matrix from the basis
of transverse-magnetic (TM) spherical vector waves into an MoM basis [63].

##  Appendix B Fundamental Bound on Realized Gain With Prescribed Input
Impedance

The upper bound on realized gain was evaluated using (6) with current vector
$\mathbf{I}$ being the solution to

$\begin{split}\underset{\mathbf{I}}{\mathrm{minimize}}\quad&-\mathbf{I}^{\mathrm{H}}\mathbf{F}^{\mathrm{H}}(\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})\mathbf{F}(\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})\mathbf{I}\\\
\mathrm{subject\,to}\quad&\mathbf{I}^{\mathrm{H}}\mathbf{Z}\mathbf{I}=\mathbf{I}^{\mathrm{H}}\mathbf{V}\\\
&\mathbf{V}^{\mathrm{H}}\mathbf{I}=\dfrac{\left|V_{\mathrm{in}}\right|^{2}}{Z_{\mathrm{in}}^{0}},\end{split}$
(10)

where the last affine constraint enforcing matching $\varGamma=0$, see (6), is
removed from the optimization using the affine transformation described in
[62]. Voltage $V_{\mathrm{in}}$, imposed on the delta-gap source together with
a real impedance $Z_{\mathrm{in}}^{0}$, is related to input power and matrix
$\mathbf{F}(\hat{\boldsymbol{d}},\hat{\boldsymbol{e}})$ gives the electric far
field.

##  Appendix C Fundamental Bound on Absorbed Power in an Uncontrollable
Subregion

This fundamental bound assumes partitioning

$\matrixquantity[\mathbf{Z}_{\mathrm{cc}}&\mathbf{Z}_{\mathrm{cu}}\\\
\mathbf{Z}_{\mathrm{uc}}&\mathbf{Z}_{\mathrm{uu}}]\matrixquantity[\mathbf{I}_{\mathrm{c}}\\\
\mathbf{I}_{\mathrm{u}}]=\matrixquantity[\mathbf{V}_{\mathrm{c}}\\\
\mathbf{V}_{\mathrm{u}}]$ (11)

into “controllable” (index c) and “uncontrollable” (index u) subregions. The
QCQP leading to the desired optimal current reads

$\displaystyle\min\limits_{\mathbf{I}_{\mathrm{c}}}$
$\displaystyle-\mathbf{I}^{\mathrm{H}}\matrixquantity[\mathbf{0}_{\mathrm{cc}}\hfil&\mathbf{0}_{\mathrm{cu}}\\\
\mathbf{0}_{\mathrm{uc}}&\mathbf{R}_{\rho,\mathrm{uu}}]\mathbf{I}$ (12)
$\displaystyle\mathrm{s.t.}$
$\displaystyle\mathbf{I}^{\mathrm{H}}\mathbf{Z}\mathbf{I}=\mathbf{I}^{\mathrm{H}}\mathbf{V}$
$\displaystyle\matrixquantity[\mathbf{Z}_{\mathrm{uc}}&\mathbf{Z}_{\mathrm{uu}}]\mathbf{I}=\mathbf{V}_{\mathrm{u}},$

where the first constraint enforces the conservation of complex power [64] and
where the last affine constraint enforcing the bottom row of the partitioned
system (11) is, as in Appendix B, removed from the optimization using the
affine transformation described in [62].

## Acknowledgement

The access to the computational infrastructure of the OP VVV funded project
CZ.02.1.01/0.0/0.0/16_019/0000765 “Research Center for Informatics” is
gratefully acknowledged.

## References

  * [1] G. Deschamps and H. Cabayan, “Antenna synthesis and solution of inverse problems by regularization methods,” _IEEE Transactions on Antennas and Propagation_ , vol. 20, no. 3, pp. 268–274, May 1972.
  * [2] M. P. Bendsoe and O. Sigmund, _Topology Optimization_ , 2nd ed. Berlin, Germany: Springer, 2004.
  * [3] M. Ohsaki, _Optimization of Finite Dimensional Structures_. CRC Press, 2011.
  * [4] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vuckovic, and A. W. Rodriguez, “Inverse design in nanophotonics,” _Nature Photonics_ , vol. 12, no. 11, pp. 659–670, 2018.
  * [5] G. Angeris, J. Vuckovic, and S. P. Boyd, “A new heuristic for physical design,” pp. 1–16, 2020. [Online]. Available: https://arxiv.org/abs/2002.05644
  * [6] J. Park and S. Boyd, “General heuristics for nonconvex quadratically constrained quadratic programming,” _arXiv preprint arXiv:1703.07870_ , 2017\.
  * [7] A. Erentok and O. Sigmund, “Topology optimization of sub-wavelength antennas,” _IEEE Trans. Antennas Propag._ , vol. 59, no. 1, pp. 58–69, 2011\.
  * [8] S. Liu, Q. Wang, and R. Gao, “MoM-based topology optimization method for planar metallic antenna design,” _Acta Mechanica Sinica_ , vol. 32, no. 6, pp. 1058–1064, Dec. 2016.
  * [9] N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund, “Giga-voxel computational morphogenesis for structural design,” _Nature_ , vol. 550, pp. 84–86, Oct. 2017.
  * [10] D. H. Wolpert and W. G. Macready, “No free lunch theorems for optimization,” _IEEE Trans. Evol. Comput._ , vol. 1, no. 1, pp. 67–82, Apr. 1997.
  * [11] M. Capek, M. Gustafsson, L. Jelinek, and P. Kadlec, “Optimal inverse design based on memetic algorithm – Part 1: Theory and implementation,” IEEE in press.
  * [12] X. Chen, C. Gu, Y. Zhang, and R. Mittra, “Analysis of partial geometry modification problems using the partitioned-inverse formula and Sherman–Morrison–Woodbury formula-based method,” _IEEE Trans. Antennas Propag_ , vol. 66, no. 10, pp. 5425–5431, October 2018.
  * [13] M. Capek, L. Jelinek, and M. Gustafsson, “Shape synthesis based on topology sensitivity,” _IEEE Trans. Antennas Propag._ , vol. 67, no. 6, pp. 3889 – 3901, June 2019.
  * [14] ——, “Inversion-free evaluation of nearest neighbors in method of moments,” _IEEE Antennas Wireless Propag. Lett._ , vol. 18, pp. 2311–2315, Apr. 2019.
  * [15] F. Jiang, S. Shen, C.-Y. Chiu, Z. Zhang, Y. Zhang, Q. S. Cheng, and R. Murch, “Pixel antenna optimization based on perturbation sensitivity analysis,” _IEEE Trans. Antennas Propag._ , vol. 70, no. 1, pp. 472–486, Jan. 2022.
  * [16] Z. Wang and S. Hum, “A new binary topology model for em surface unit cell design optimization,” in _2022 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting (AP-S/URSI)_ , 2022, pp. 67–68.
  * [17] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, _Introduction to Algorithms_ , 3rd ed. Massachusetts Institute of Technology, 2009.
  * [18] D. Simon, _Evolutionary Optimization Algorithms_. John Wiley & Sons, 2013.
  * [19] G. C. Onwubolu and B. V. Babu, _New Optimization Techniques in Engineering_. Berlin, Germany: Springer, 2004.
  * [20] O. Sigmund, “On the usefulness of non-gradient approaches in topology optimization,” _Struct. Multidisc. Optim._ , vol. 43, pp. 589–596, 2011\.
  * [21] S. Katoch, S. S. Chauhan, and V. Kumar, “A review on genetic algorithm: Past, present, and future,” _Multimedia Tools and Applications_ , vol. 80, no. 5, pp. 8091–8126, Oct. 2020.
  * [22] Local Optimization Based on Topology Sensitivity. [Online]. Available: https://github.com/MiloslavCapek/TopologySensitivity
  * [23] (2018) The Matlab. The MathWorks. [Online]. Available: www.mathworks.com
  * [24] (2019) Antenna Toolbox for MATLAB (AToM). Czech Technical University in Prague. www.antennatoolbox.com. [Online]. Available: {www.antennatoolbox.com}
  * [25] (2017) Fast Optimization Procedures (FOPS). Czech Technical University in Prague. [Online]. Available: www.antennatoolbox.com/fops
  * [26] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” _IEEE Trans. Antennas Propag._ , vol. 30, no. 3, pp. 409–418, May 1982.
  * [27] J. A. De Loera, J. Rambau, and F. Santos, _Triangulations – Structures for Algorithms and Applications_. Berlin, Germany: Springer, 2010.
  * [28] D. A. Dunavant, “High degree efficient symmetrical gaussian quadrature rules for the triangle,” _International Journal for Numerical Methods in Engineering_ , vol. 21, pp. 1129–1148, 1985.
  * [29] (2021) RCI cluster. Research Center for Informatics, CTU FEE, Czech Republic. [Online]. Available: https://login.rci.cvut.cz/wiki/start
  * [30] K. Schab, L. Jelinek, M. Capek, C. Ehrenborg, D. Tayli, G. A. E. Vandenbosch, and M. Gustafsson, “Energy stored by radiating systems,” _IEEE Access_ , vol. 6, pp. 10 553–10 568, 2018.
  * [31] M. Capek, L. Jelinek, K. Schab, M. Gustafsson, B. L. G. Jonsson, F. Ferrero, and C. Ehrenborg, “Optimal planar electric dipole antennas: Searching for antennas reaching the fundamental bounds on selected metrics,” _IEEE Antennas and Propagation Magazine_ , vol. 61, no. 4, pp. 19–29, Aug. 2019.
  * [32] A. D. Yaghjian and S. R. Best, “Impedance, bandwidth and Q of antennas,” _IEEE Trans. Antennas Propag._ , vol. 53, no. 4, pp. 1298–1324, Apr. 2005\.
  * [33] J. L. Volakis, C. Chen, and K. Fujimoto, _Small Antennas: Miniaturization Techniques & Applications_. McGraw-Hill, 2010.
  * [34] M. Capek, M. Gustafsson, and K. Schab, “Minimization of antenna quality factor,” _IEEE Trans. Antennas Propag._ , vol. 65, no. 8, pp. 4115–4123, Aug, 2017.
  * [35] S. R. Best, “Electrically small resonant planar antennas,” _IEEE Antennas Propag. Mag._ , vol. 57, no. 3, pp. 38–47, June 2015.
  * [36] K. Fujimoto and H. Morishita, _Modern Small Antennas_. Cambridge, Great Britain: Cambridge University Press, 2013\.
  * [37] S. R. Best, “Low Q electrically small linear and elliptical polarized spherical dipole antennas,” _IEEE Trans. Antennas Propag._ , vol. 53, no. 3, pp. 1047–1053, 2005.
  * [38] M. Gustafsson, D. Tayli, C. Ehrenborg, M. Cismasu, and S. Norbedo, “Antenna current optimization using MATLAB and CVX,” _FERMAT_ , vol. 15, no. 5, pp. 1–29, May–June 2016.
  * [39] M. Capek and L. Jelinek, “Optimal composition of modal currents for minimal quality factor Q,” _IEEE Trans. Antennas Propag._ , vol. 64, no. 12, pp. 5230–5242, Dec. 2016.
  * [40] L. Jelinek and M. Capek, “Optimal currents on arbitrarily shaped surfaces,” _IEEE Trans. Antennas Propag._ , vol. 65, no. 1, pp. 329–341, Jan. 2017.
  * [41] M. Capek, V. Neuman, J. Tucek, L. Jelinek, and M. Gustafsson, “Topology optimization of electrically small antennas with shape regularity constraints,” in _Proceedings of the 15th European Conference on Antennas and Propagation (EUCAP)_ , 2021.
  * [42] K. Svanberg, “The method of moving asymptotes—a new method for structural optimization,” _International Journal for Numerical Methods in Engineering_ , vol. 24, no. 2, pp. 359–373, Feb. 1987.
  * [43] T. E. Bruns and D. A. Tortorelli, “Topology optimization of non-linear elastic structures and compliant mechanisms,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 190, no. 26-27, pp. 3443–3459, Mar. 2001.
  * [44] F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” _Structural and Multidisciplinary Optimization_ , vol. 43, no. 6, pp. 767–784, Dec. 2010.
  * [45] B. Yang, J. Zhou, and J. J. Adams, “A Shape-First, Feed-Next Design Approach For Compact Planar MIMO Antennas,” _Progress In Electromagnetics Research M_ , vol. 77, pp. 157–165, 2019.
  * [46] S. R. Best, “The radiation properties of electrically small folded spherical helix antennas,” _IEEE Trans. Antennas Propag._ , vol. 52, no. 4, pp. 953–960, 2004.
  * [47] H. L. Thal, “New radiation Q limits for spherical wire antennas,” _IEEE Trans. Antennas Propag._ , vol. 54, no. 10, pp. 2757–2763, Oct. 2006\.
  * [48] ——, “Gain and Q bounds for coupled TM-TE modes,” _IEEE Trans. Antennas Propag._ , vol. 57, no. 7, pp. 1879–1885, July 2009.
  * [49] M. Gustafsson, M. Capek, and K. Schab, “Tradeoff between antenna efficiency and Q-factor,” _IEEE Trans. Antennas Propag._ , vol. 67, no. 4, pp. 2482–2493, April 2019.
  * [50] D. M. Pozar, _Microwave Enginnering_ , 4th ed. Wiley, 2011.
  * [51] M. Ehrgott, _Multicriteria Optimization_. Springer.
  * [52] L. J. Cohon, _Multiobjective Programming and Planning_ , 1st ed. Academic Press, 1978.
  * [53] C. A. Balanis, _Modern Antenna Handbook_. Wiley, 2008.
  * [54] J. D. Jackson, _Classical Electrodynamics_ , 3rd ed. Wiley, 1998.
  * [55] C. A. Balanis, _Antenna Theory Analysis and Design_ , 3rd ed. Wiley, 2005.
  * [56] Y. Rahmat-Samii and E. Michielssen, Eds., _Electromagnetic Optimization by Genetic Algorithms_. Wiley, 1999.
  * [57] G. Strang, _Introduction to Applied Mathematics_. Cambridge University Press, 1986.
  * [58] G. Bendall and F. Margot, “Greedy-type resistance of combinatorial problems,” _Discrete Optimization_ , vol. 3, no. 4, pp. 288–298, Dec. 2006.
  * [59] R. L. Haupt and D. H. Werner, _Genetic Algorithms in Electromagnetics_. John Wiley & Sons, 2007\.
  * [60] K. Deb, _Multi-Objective Optimization using Evolutionary Algorithms_. New York, United States: Wiley, 2001.
  * [61] M. Cismasu and M. Gustafsson, “Antenna bandwidth optimization with single frequency simulation,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 3, pp. 1304–1311, 2014.
  * [62] J. Liska, L. Jelinek, and M. Capek, “Fundamental Bounds to Time-Harmonic Quadratic Metrics in Electromagnetism: Overview and Implementation,” 2021.
  * [63] V. Losenicky, L. Jelinek, M. Capek, and M. Gustafsson, “Method of moments and T-matrix hybrid,” _IEEE Transactions on Antennas and Propagation_ , vol. 70, no. 5, pp. 3560–3574, 2022.
  * [64] M. Gustafsson, K. Schab, L. Jelinek, and M. Capek, “Upper bounds on absorption and scattering,” _New Journal of Physics_ , vol. 22, no. 7, p. 073013, sep 2020.

|  Miloslav Capek (M’14, SM’17) received the M.Sc. degree in Electrical
Engineering 2009, the Ph.D. degree in 2014, and was appointed Associate
Professor in 2017, all from the Czech Technical University in Prague, Czech
Republic. He leads the development of the AToM (Antenna Toolbox for Matlab)
package. His research interests are in the area of electromagnetic theory,
electrically small antennas, numerical techniques, fractal geometry, and
optimization. He authored or co-authored over 120 journal and conference
papers. Dr. Capek is Associate Editor of IET Microwaves, Antennas &
Propagation. He was a regional delegate of EurAAP between 2015 and 2020. He
received the IEEE Antennas and Propagation Edward E. Altshuler Prize Paper
Award 2023.
---|---
|  Lukas Jelinek received his Ph.D. degree from the Czech Technical University
in Prague, Czech Republic, in 2006. In 2015 he was appointed Associate
Professor at the Department of Electromagnetic Field at the same university.
His research interests include wave propagation in complex media,
electromagnetic field theory, metamaterials, numerical techniques, and
optimization.
---|---
|  Petr Kadlec (M’13) received the Ph.D. degree in electrical engineering from
the Brno University of Technology (BUT), Brno, Czech Republic, in 2012. He is
currently an Associate Professor with the Department of Radioelectronics, BUT.
His research interests include global optimization methods and computational
methods in electromagnetics. He is a leading developer of the FOPS (Fast
Optimization ProcedureS) MATLAB software package.
---|---
|  Mats Gustafsson received the M.Sc. degree in Engineering Physics 1994, the
Ph.D. degree in Electromagnetic Theory 2000, was appointed Docent 2005, and
Professor of Electromagnetic Theory 2011, all from Lund University, Sweden. He
co-founded the company Phase holographic imaging AB in 2004. His research
interests are in scattering and antenna theory and inverse scattering and
imaging. He has written over 100 peer reviewed journal papers and over 100
conference papers. Prof. Gustafsson received the IEEE Schelkunoff Transactions
Prize Paper Award 2010, the IEEE Uslenghi Letters Prize Paper Award 2019, and
best paper awards at EuCAP 2007 and 2013. He served as an IEEE AP-S
Distinguished Lecturer for 2013-15.
---|---
  *[DOF]: degrees-of-freedom
  *[MoM]: method of moments
  *[PEC]: perfect electric conductor
  *[RWG]: Rao-Wilton-Glisson
  *[QCQP]: quadratically constrained quadratic program


version of the 2D diffeomorphism symmetry transformations. Using the
theoretical potential of MBTSA, we have established that the CF-type
restrictions, in the case of the 2D diffeomorphism invariant theory, is the
generalization of the 1D version of the CF-type restriction that exists for
the 1D (non-)interacting systems of the (non-)SUSY and (non-)relativistic
particles which have been discussed in our present endeavor. It is pertinent
to point out that the 1D and 2D version of the CF-type restriction(s) are the
special cases of the D-dimensional CF-type restrictions that exist for the
D-dimensional diffeomorphism invariant theory [7, 8]. This establishes the
universal nature of the CF-type restriction.

We dwell a bit now on the existence of CF-type restriction(s) in the context
of BRST formalism and its importance therein. As we have seen in our present
endeavor that the gauge symmetries are generated by the first-class
constraints in the terminology of the Dirac prescription scheme for the
classification of constraints [13, 14]. These gauge symmetries can be
exploited within the framework of BRST formalism for the covariant canonical
quantization of the gauge theories. Thus, the most modern and fundamental
definition of a classical gauge theory is the observation that a theory that
is endowed with the first-class constraints (which generate the local gauge
symmetry transformations for the gauge theory under consideration). In exactly
similar fashion, we state that when a gauge and/or diffeomorphism invariant
theory is quantized (within the framework of BRST formalism), the
characteristic feature of the BRST-quantized theory is the observation that
(non-) trivial CF-type restriction exists on the quantized gauge and/or
diffeomorphism invariant theory. For the gauge theory, discussed within the
framework of BRST formalism, the CF-type restriction(s) are not universal and
their form and number might change from a specific system to the other systems
(see, e.g. [19, 20, 21, 23] for details). However, for the D-dimensional
diffeomorphism invariant theory, we have been able to show that the form of
the CF-type restrictions is:
$B_{\mu}+\bar{B}_{\mu}+i\,(\bar{C}^{\rho}\,\partial_{\rho}\,C_{\mu}+C^{\rho}\,\partial_{\rho}\,\bar{C}_{\mu})=0$
(with $\mu,\,\rho=0,1,2,...,D-1$). Here the bosonic fields
$(B_{\mu},\bar{B}_{\mu})$ are the Nakanishi-Lautrup type auxiliary fields and
the fermionic (anti-)ghost fields are denoted by $(\bar{C}_{\mu})C_{\mu}$. It
is straightforward to note that, in our present endeavor, we have obtained the
1D and 2D versions of this (anti-)BRST invariant universal CF-type
restrictions (valid for the D-dimensional diffeomorphism invariant theory) as
we have focused on the 1D and 2D diffeomorphism invariant theories.

One of the hallmarks of the gravitational as well as (super)string theories is
the diffeomorphism invariance. Thus, our present article is a precursor to the
BRST and related superfield approach to the diffeomorphism invariant models of
physical interest in the physical four (3+1)-dimensions of spacetime. This is
the future perspective and scope of our present review-cum-research article.
We have plans to pursue this direction in our future endeavors. The results of
our research works will be reported elsewhere [40].

Acknowledgments

One of us (AKR) gratefully acknowledges the financial support from the
Institution of Eminence (IoE) Research Incentive Grant of PFMS Scheme No.
3254-World Class Institutions to which Banaras Hindu University, Varanasi,
belongs. The present investigation has been carried out under the above
Research Incentive Grant which has been launched by the Govt. of India to
promote academic excellence in research.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Appendix A: On the Singular Nature of the Lagrangian for the Non-SUSY 1D
Diffeomorphism Invariant Scalar Relativistic Particle

To corroborate our assertion that our starting Lagrangian
$L_{0}=m\,\sqrt{{\dot{x}}^{\mu}\,\dot{x}_{\mu}}\equiv\sqrt{\dot{x}_{0}^{2}-\dot{x}_{i}^{2}}\;(i=1,2,...D-1)$
for the free scalar relativistic particle is singular in any arbitrary
D-dimension of target spacetime, we take up the three (2 + 1)-dimensional (3D)
flat Minkowski spacetime manifold as the target space for which we have:
$\alpha,\beta,\gamma...\mu,\nu,\lambda...=0,1,2$ and $i,j,k...=1,2$. In this
case, we have the following $3\times 3$ matrix (see, e.g. [15, 14] for
details)

$\displaystyle
W_{\alpha\beta}\;=\;\large{\begin{pmatrix}W_{00}\;&\;W_{01}\;&\;W_{02}\\\
W_{10}\;&\;W_{11}\;&\;W_{12}\\\
W_{20}\;&\;W_{21}\;&\;W_{22}\end{pmatrix}\quad\equiv\quad\begin{pmatrix}\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{0}\,\partial\,\dot{x}^{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{0}\,\partial\,\dot{x}^{1}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{0}\,\partial\,\dot{x}^{2}}}\\\
\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{1}\,\partial\,\dot{x}^{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{1}\,\partial\,\dot{x}^{1}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{1}\,\partial\,\dot{x}^{2}}}\\\
\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{2}\,\partial\,\dot{x}^{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{2}\,\partial\,\dot{x}^{1}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{2}\,\partial\,\dot{x}^{2}}}\end{pmatrix},~{}~{}~{}~{}~{}~{}~{}~{}(A.1)}$

where $x^{\alpha}=(x^{0},x^{1},x^{2})\equiv(x_{0},\underline{x},y)$ and
$x_{\alpha}=(x_{0},x_{1},x_{2})\equiv(x_{0},-\,\underline{x},-\,y)$ so that we
have: $x_{\mu}\,x^{\mu}=x\cdot x=x^{2}=x_{0}^{2}-x_{1}^{2}-x_{2}^{2}\equiv
x_{0}^{2}-\underline{x}^{2}-y^{2}$. The sign difference in $x^{\mu}$ and
$x_{\mu}$ is due to our choice of the 3D metric $\eta_{\mu\nu}$ = diag
$(+1,-1,-1)\equiv\eta^{\mu\nu}$ and $x^{1}=\underline{x}=-\,x_{1}$ has been
taken into account to differentiate between
$x^{2}=x_{0}^{2}-\underline{x}^{2}-y^{2}$ and
$x_{1}^{2}\equiv\underline{x}^{2}$. As a result of this choice of notations,
we have the explicit form of the above $3\times 3$ matrix as:

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{\alpha\beta}\;=\;\large{\begin{pmatrix}\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}_{0}\,\partial\,\dot{x}_{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}_{0}\,\partial\,\dot{{\underline{x}}}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}_{0}\,\partial\,\dot{y}}}\\\
\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{{\underline{x}}}\,\partial\,\dot{x}_{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{{\underline{x}}}\,\partial\,\dot{{\underline{x}}}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{{\underline{x}}}\,\partial\,\dot{y}}}\\\
\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{y}\,\partial\,\dot{x}_{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{y}\,\partial\,\dot{{\underline{x}}}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{y}\,\partial\,\dot{y}}}\end{pmatrix}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(A.2)}$

The substitutions of all the entries into the above $3\times 3$ matrix (with
proper considerations of the precise notations) lead to the following

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{\alpha\beta}\;=\;\frac{m}{(\dot{x}^{2})^{3/2}}\large{\begin{pmatrix}-\,(\dot{{\underline{x}}^{2}}+\dot{y}^{2})\;&\;\dot{x}_{0}\,\dot{{\underline{x}}}\;&\;\dot{x}_{0}\,\dot{y}\\\
\dot{x}_{0}\,\dot{{\underline{x}}}\;&\;(\dot{y}^{2}-\dot{x}_{0}^{2})\;&\;-\,\dot{{\underline{x}}}\,\dot{y}\\\
\dot{x}_{0}\,\dot{y}\;&\;-\,\dot{{\underline{x}}}\,\dot{y}\;&\;(\dot{{\underline{x}}^{2}}-\dot{x}_{0}^{2})\end{pmatrix},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(A.3)}$

where we have taken into account:
$\dot{x}^{2}-{\dot{x}}_{0}^{2}=-\,(\underline{\dot{x}}^{2}+\dot{y}^{2}).$ Now,
at this stage, it is straightforward to note that we have the following
expression for the determinant, namely;

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mbox{det}\;W_{\alpha\beta}=[m^{3}/(\dot{x}^{2})^{9/2}]\,[\dot{x}_{0}^{2}\,\\{(\dot{{\underline{x}}^{2}}+\dot{y}^{2})-(\dot{{\underline{x}}^{2}}+\dot{y}^{2})\\}]=0.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(A.4)$

We end this Appendix with the final remark that we have shown that det
$W_{\alpha\beta}=0$ for the $(2\times 2)$ and ($3\times 3$) matrices, in our
present endeavor, for the $D=2,3$ dimensional Minkowskian flat spacetime
manifolds that are taken as the target space, respectively. However, this
observation is true for any arbitrary dimension of the flat Minkowskian
spacetime manifold that can be chosen for the target space.

Appendix B: On the Singular Nature of the Starting Lagrangian for the Free
Massive Spinning Relativistic Particle

We take into consideration the three (2 + 1)-dimensional (3D) flat Minkowskian
spacetime manifold as the target space to show that, even in the 3D case, the
determinant of the $3\times 3$ matrix $W_{\alpha\beta}$ turns out to be zero.
For this purpose, we focus on the starting Lagrangian $L_{0}$ of Eq. (31) and
compute the following $3\times 3$ matrix (see, e.g. [15, 14] for details)

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}W_{\alpha\beta}\;=\;\large{\begin{pmatrix}\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{0}\,\partial\,\dot{x}^{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{0}\,\partial\,\dot{x}^{1}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{0}\,\partial\,\dot{x}^{2}}}\\\
\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{1}\,\partial\,\dot{x}^{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{1}\,\partial\,\dot{x}^{1}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{1}\,\partial\,\dot{x}^{2}}}\\\
\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{2}\,\partial\,\dot{x}^{0}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{2}\,\partial\,\dot{x}^{1}}}\;&\;\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{2}\,\partial\,\dot{x}^{2}}}\end{pmatrix}}\quad\equiv\quad\frac{{\partial^{2}\,L_{0}}}{{\partial\,\dot{x}^{\alpha}\,\partial\,\dot{x}^{\beta}}},~{}~{}~{}~{}~{}~{}~{}~{}~{}(B.1)$

where we have to take the help of explicit form of $W_{\alpha\beta}$ of Eq.
(33) for computations.

At this stage, due to our choice of the metric tensor
$\eta_{\mu\nu}=(+1,-\,1,\,-\,1)\equiv\eta^{\mu\nu}$, we have the explicit form
of the notations:
$\dot{X}^{\mu}=(\dot{X}^{0},\,\dot{X}^{1},\,\dot{X}^{2})=(\dot{X}_{0},\,\dot{\underline{X}},\,\dot{Y})$
and
$\dot{X}_{\mu}=(\dot{X}_{0},\,\dot{X}_{1},\,\dot{X}_{2})=(\dot{X}_{0},\,-\,\dot{\underline{X}},\,-\,\dot{Y})$
so that the symbol $A=\dot{X}^{\mu}\,\dot{X}_{\mu}$ is expanded as:
$A=\dot{X}^{2}\equiv\dot{X}_{0}^{2}-\dot{X}_{1}^{2}-\dot{X}_{2}^{2}=\dot{X}_{0}^{2}-\dot{\underline{X}}^{2}-\dot{Y}^{2}$
which lead to the following $3\times 3$ matrix, namely;

$\displaystyle
W_{\alpha\beta}\;=\;\large{\begin{pmatrix}W_{00}\;&\;W_{01}\;&\;W_{02}\\\
W_{10}\;&\;W_{11}\;&\;W_{12}\\\ W_{20}\;&\;W_{21}\;&\;W_{22}\end{pmatrix}}$
$\displaystyle\quad~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\equiv\quad\frac{m}{(A)^{3/2}}{\begin{pmatrix}-(\dot{\underline{X}}^{2}+\dot{Y}^{2})\;&\;\dot{X}_{0}\,\dot{\underline{X}}\;&\;\dot{X}_{0}\,\dot{Y}\\\
\dot{X}_{0}\,\dot{\underline{X}}\;&\;(\dot{Y}^{2}-\dot{X}_{0}^{2})\;&\;-\,\dot{\underline{X}}\,\dot{Y}\\\
\dot{X}_{0}\,\dot{Y}\;&\;-\,\dot{\underline{X}}\,\dot{Y}\;&\;(\dot{\underline{X}}^{2}-\dot{X}_{0}^{2})\end{pmatrix}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(B.2)}$

It should be noted that we have chosen the notations of $X^{\mu}$ and
$X_{\mu}$ in the component forms such that there is a difference between
$A=\dot{X}^{\mu}\,\dot{X}_{\mu}=\dot{X}^{2}$ and
${\dot{\underline{X}}}^{2}={\dot{X}}_{1}\,\dot{X}_{1}$. In other words, we
have:
$\dot{X}^{2}=\dot{X}_{0}^{2}-\dot{X}_{1}^{2}-\dot{X}_{2}^{2}=\dot{X}_{0}^{2}-\dot{\underline{X}}^{2}-\dot{Y}^{2}$
and $\dot{X}_{1}^{2}={\dot{\underline{X}}}^{2}$, respectively. The explicit
computations of the above determinant leads to the following:

$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mbox{det}\,W_{\alpha\beta}=\frac{m^{3}}{(A)^{9/2}}\,\Big{[}\dot{X}_{0}^{2}\;\\{(\dot{\underline{X}}^{2}+\dot{Y}^{2})-(\dot{\underline{X}}^{2}+\dot{Y}^{2})\\}\Big{]}=0.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(B.3)$

Thus, the starting Lagrangian $(L_{0})$ of (31) is a singular Lagrangian. We
have proven the singular nature of $(L_{0})$ in two and three dimensional flat
Minkowskian target spacetime manifold. However, our statement is true for any
arbitrary D-dimensional flat Minkowskian target spacetime manifold. It is also
pretty obvious that $(L_{0})$ is singular w.r.t. the fermionic ($\chi^{2}=0$)
variable $\chi$ because it has no generalized velocity ($\dot{\chi}$) in it.
In fact, in the whole theory (that is described by $L_{0}$), we have the
existence of two constraints: $p^{2}-m^{2}\approx 0$ and
$p_{\mu}\,\psi^{\mu}-m\,\psi_{5}\approx 0$ which can be explicitly checked.

Appendix C: On the Derivation of the (Anti-)BRST Symmetries for Other
Variables of Our Interacting 1D Diffeomorphism Invariant Theory

In this Appendix, we focus on the derivation of the (anti-)BRST symmetries of
the rest of the variables of the 1D diffeomorphism invariant interacting
scalar relativistic particle besides the ones that have been derived by using
the theoretical strength of the MBTSA to BRST formalism. We exploit the
requirements of (i) the off-shell nilpotency, (ii) the absolute
anticommutativity, and (iii) the invariance of the CF-type restriction under
the (anti-)BRST symmetry transformations for this purpose. In this regard, it
is pertinent to point out that the application of MBTSA to BRST formalism
leads to the derivation of the following (anti-)BRST symmetry transformations
and CF-type restriction

$\displaystyle s_{b}x_{\mu}=C\,\dot{x}_{\mu},\qquad
s_{b}p_{\mu}=C\,\dot{p}_{\mu},\qquad s_{ab}x_{\mu}=\bar{C}\,\dot{x}_{\mu},$
$\displaystyle s_{ab}p_{\mu}=\bar{C}\,\dot{p}_{\mu},\qquad
B+\bar{B}+i\,(\bar{C}\,\dot{C}-\dot{\bar{C}}\,C)=0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(C.1)$

where we have also used the standard (anti-)BRST symmetry transformations:
$s_{b}\bar{C}=i\,B$ and $s_{ab}C=i\,\bar{B}$. It is straightforward to note
that the requirement of the off-shell nilpotency of all the above (anti-)BRST
symmetry transformations [i.e. $s_{(a)b}^{2}=0$] leads to:

$s_{b}B=0,\qquad s_{ab}\bar{B}=0,\qquad s_{b}C=C\,\dot{C},\qquad
s_{ab}\bar{C}=\bar{C}\,\dot{\bar{C}}.$ $None$

On the other hand, the requirement of the (anti-)BRST invariance of the CF-
type restriction: [i.e. $B+\bar{B}+i\,(\bar{C}\,\dot{C}-\dot{\bar{C}}\,C)=0$]
leads to the following

$s_{b}\,\bar{B}=\dot{\bar{B}}\,C-\bar{B}\,\dot{C},\qquad
s_{ab}\,B=\dot{B}\,\bar{C}-B\,\dot{\bar{C}},$ $None$

provided we use the validity of the CF-type restriction on our theory. It is
straightforward to check that the (anti-)BRST symmetry transformations, that
have been derived in (C.2) and (C.3), are off-shell nilpotent
$[s_{(a)b}^{2}=0]$ of order two. Furthermore, it is interesting to point out
that the requirement of the following absolute anticommutativity property on
the (anti-)ghost variables, namely;

$\\{s_{b},\,s_{ab}\\}\,\bar{C}=0,\qquad\qquad\\{s_{b},\,s_{ab}\\}\,C=0,$
$None$

also leads to the derivation of the (anti-)BRST symmetry symmetry
transformations for the Nakanishi-Lautrup variables $(\bar{B},\,B)$ quoted in
(C.3). In addition, it can be checked that whatever (anti-)BRST symmetry
transformations are derived, so far, are proper in the sense that they satisfy
(i) the off-shell nilpotency property, and (ii) the absolute anticommutativity
property. The invariance of the CF-type restriction, under the above
(anti-)BRST symmetry transformations, is automatic and natural which
demonstrates the physical nature of this restriction (i.e. CF-type
restriction) on our theory.

We end this Appendix with the remark that the theoretical strength of the
MBTSA to BRST formalism plus the requirements of the off-shell nilpotency,
absolute anticommutativity and invariance of the CF-type restrictions are not
good enough to lead to the derivation of the (anti-)BRST symmetry
transformations for the einbein variable $e\,(\tau)$ for the
reparameterization invariant model of an interacting scalar relativistic
particle with the electromagnetic field in the background. However, we know
that the einbein variable $e(\tau)$ is like a “gauge” variable under the
classical gauge symmetry transformations and transforms like a “scalar
density” under the classical diffeomorphism transformations. Hence, the
following generalizations of the classical reparameterization symmetries to
their counterparts quantum (anti-)BRST symmetries, namely;

$\delta_{r}\,e=\frac{d}{d\,\tau}\,(\varepsilon\,e)\;\;\longrightarrow\;\;s_{b}\,e=\frac{d}{d\,\tau}\,(C\,e),$
$\delta_{r}\,e=\frac{d}{d\,\tau}\,(\varepsilon\,e)\;\;\longrightarrow\;\;s_{ab}\,e=\frac{d}{d\,\tau}\,(\bar{C}\,e),$
$None$

satisfy all the requirements of the BRST formalism in the sense that the off-
shell nilpotency and absolute anticommutativity properties are respected by
the (anti-)BRST symmetry transformations quoted in (C.5) provided we take into
account the sanctity and validity of the CF-type restriction [i.e.
$B+\bar{B}+i\,(\bar{C}\,\dot{C}+\dot{\bar{C}}\,C)=0$].

## References

  * [1] M. B. Green, J. H. Schwarz, E. Witten, Superstring Theory, Vols. 1 and 2,
Cambridge University Press, Cambridge (1987)

  * [2] J. Polchinski, String Theory, Vols. 1 and 2,
Cambridge University Press, Cambridge (1998)

  * [3] D. Lust, S. Theisen, Lectures in String Theory, Springer-Verlag, New York (1989)
  * [4] K. Becker, M. Becker, J. H. Schwarz, String Theory and M-Theory,
Cambridge University Press, Cambridge (2007)

  * [5] D. Rickles, A Brief History of String Theory From Dual Models to M-Theory,
Springer, Germany (2014)

  * [6] M. Kato, K. Ogawa, Covariant Quantization of String Based on BRS Invariance,
Nucl. Phys.  B 212, 443 (1983)

  * [7] L. Bonora, BRST and Supermanifolds, Nucl. Phys. B 912, 103 (2016)
  * [8] L. Bonora, R. P. Malik, BRST and Superfield Formalism–A Review,
Universe 7, 280 (2021)

  * [9] L. Bonora, M. Tonin, Superfield Formulation of Extended BRS Symmetry,
Phys. Lett. B 98, 48 (1981)

  * [10] L. Bonora, P. Pasti, M. Tonin, Extended BRS Symmetry in Non-Abelian Gauge
Theories, Nuovo Cimento A 64, 307 (1981)

  * [11] L. Bonora, P. Pasti, M. Tonin, Superspace Approach to Quantum Gauge Theories, Ann. Phys. 144, 15 (1982)
  * [12] G. Curci, R. Ferrari, Slavnov Transformations and Supersymmetry,
Phys. Lett. B 63, 91 (1976)

  * [13] P. A. M. Dirac, Lectures on Quantum Mechanics (Belfer Graduate School of Science), Yeshiva University Press, New York (1964)
  * [14] K. Sundermeyer, Constraint Dynamics, Lecture Notes in Physics,
Vol. 169, Springer-Verlag, Berlin (1982)

  * [15] E. C. G. Sudarshan, N. Mukunda, Classical Dynamics: A Modern Perspective,
Wiley, New York (1972)

  * [16] P. Mitra, R. Rajaraman, New Results on Systems with Second-Class Constraints,
Ann. Phys. 203, 137 (1990)

  * [17] P. Mitra, R. Rajaraman, Gauge-Invariant Reformulation of Theories with Second-Class Constraints, Ann. Phys. 203, 157 (1990)
  * [18] D. M. Gitman, I. V. Tyutin, Quantization of Fields with Constraints,
Springer-Verlag, Berlin Heidelberg (1990)

  * [19] L. Bonora, R. P. Malik, BRST, Anti-BRST and Gerbes,
Phys. Lett. B 655, 75 (2007)

  * [20] L. Bonora, R. P. Malik, BRST, Anti-BRST and Their Geometry,
J. Phys. A: Math. Theor. 43, 375403 (2010)

  * [21] A. K. Rao, A. Tripathi, R. P. Malik, Supervariable and BRST Approaches to a Reparameterization Invariant Nonrelativistic System,
Adv. High Energy Phys. 2021, 5593434 (2021)

  * [22] B. Chauhan, A. K. Rao, A. Tripathi, R. P. Malik, 1D Diffeomorphism Invariant Model of a Free Scalar Relativistic Particle: Supervariable and BRST Approaches,
Int. J. Mod. Phys. 37, 2250164 (2022)

  * [23] A. Tripathi, B. Chauhan, A. K. Rao, R. P. Malik, Reparameterization Invariant Model of a Supersymmetric System: BRST and Supervariable Approaches,
Adv. High Energy Phys. 2021, 2056629 (2021)

  * [24] R. P. Malik, Nilpotent Symmetries of a Model of 2D Diffeomorphism Invariant Theory: BRST Approach, Adv. High Energy Phys. 2022, 8155214 (2022)
  * [25] A. Tripathi, A. K. Rao, R. P. Malik, Superfield Approaches to a Model of Bosonic String: Curci-Ferrari-Type Restrictions, Adv. High Energy Phys. 2022, 9505924 (2022)
  * [26] R. Brunetti, K. Fredenhagen, Time of Occurrence Observable in Quantum Mechanics, Phys. Rev. A 66, 044101 (2002)
  * [27] R. Banerjee, B. Chakraborty, S. Gangopadhyay, Noncommutativity and Reparametrization Symmetry, J. Phys. A: Math. Gen. 38, 957 (2005)
  * [28] R. P. Malik, Noncommutativity in a Simple Toy Model,
Int. J. Mod. Phys. A 22, 2660 (2007)

  * [29] D. Nemschansky, C. Preitschopf, M. Weinstein, A BRST Primer,
Ann. Phys. 183, 226 (1988)

  * [30] L. Brink, P. Di Vecchia, P. Howe, A Lagrangian Formulation of the Classical and Quantum Dynamics of Spinning Particles, Nucl. Phys. B 118, 76 (1977)
  * [31] L. Brink, S. Deser, B. Zumino, P. Di Vecchia, P. Howe, Local Supersymmetry for Spinning Particles, Phys. Lett. B 64, 435 (1976)
  * [32] A. Shukla, S. Krishna, R. P. Malik, Supersymmetrization of Horizontality Condition: Nilpotent Symmetries for a Free Spinning Relativistic Particle,
Eur. Phys. J. C 72, 2188 (2012)

  * [33] R. P. Malik, Interacting Relativistic Particle: Time-Space Noncommutativity and Symmetries, Int. J. Mod. Phys. A 22, 1053 (2007)
  * [34] A. Shukla, N. Srinivas, R. P. Malik, Nilpotent Symmetries of a 4D Model of Hodge Theory: Augmented (Anti-)Chiral Superfield Formalism, Ann. Phys. 394, 98 (2018)
  * [35] B. Chauhan, S. Kumar, A. Tripathi, R. P. Malik, Modified 2D Proca Theory: Revisited Under BRST and (Anti-)Chiral Superfield Formalisms,
Adv. High Energy Phys. 2020, 3495168 (2020)

  * [36] S. Kumar, B. Chauhan, A. Tripathi, R. P. Malik, Massive 4D Abelian 2-Form Theory: Nilpotent Symmetries from the (Anti-)Chiral Superfield Approach,
Int. J. Mod. Phys. 37, 2250003 (2022)

  * [37] A. M. Polyakov, Quantum Geometry of Bosonic Strings,
Phys. Lett. B 103, 207 (1981)

  * [38] A. M. Polyakov, Quantum Geometry of Fermionic Strings,
Phys. Lett. B 103, 211 (1981)

  * [39] T. Kugo, S. Uehara, General Procedure of Gauge Fixing Based on BRS Invariance Principle, Nucl. Phys. B 197, 378 (1982)
  * [40] R. P. Malik, etal., in preparation