# Density-Based Topology Optimization in
Method of Moments: Q-factor Minimization

Jonas Tucek, Miloslav Capek, Lukas Jelinek, and Ole Sigmund Manuscript
received XXX; revised XXX. This work was supported by the Czech Science
Foundation under project No. 21-19025M, by the Czech Technical University in
Prague under project SGS22/162/OHK3/3T/13, and by financial support from the
Villum Foundation through the Villum Investigator Project InnoTop. J. Tucek,
M. Capek and L. Jelinek are with the Czech Technical University in Prague,
Prague, Czech Republic (e-mails: jonas.tucek; miloslav.capek;
lukas.jelinek@fel.cvut.cz).O. Sigmund is with Technical University of Denmark,
Lyngby, Denmark (e-mail: olsi@dtu.dk).

###### Abstract

Classical gradient-based density topology optimization is adapted for method-
of-moments numerical modeling to design a conductor-based system attaining the
minimal antenna Q-factor evaluated via an energy stored operator. Standard
topology optimization features are discussed, e.g., interpolation scheme and
density and projection filtering. The performance of the proposed technique is
demonstrated in a few examples in terms of the realized Q-factor values and
necessary computational time to obtain a design. The optimized designs are
compared to the fundamental bound and well-known empirical structures. The
presented framework can provide a completely novel design, as presented in the
second example.

###### Index Terms:

Antennas, Topology optimization, numerical methods, Q-factor.

##  I Introduction

Wireless communication systems impose more and more requirements on the
performance of antenna systems, e.g., compactness, efficiency, or large
operational bandwidth. Standard antenna design methodologies [1] are often
unsuccessful in delivering an antenna satisfying these demands which are
contradictory in nature [2], [3], and designers are forced to utilize advanced
procedures based on numerical modeling.

Antenna synthesis attempts to determine the spatial material distribution and
feeding network to meet performance requirements and satisfy constraints. This
paper focuses solely on optimizing material distribution with predetermined
and fixed feeding properties to reduce the complexity of the optimization
problem.

Antenna engineers often exploit their previous empirical knowledge when fully
describing a structure with a set of parameters which are further optimized by
simple parametric sweeps or surrogate-based modeling [4] in combination with
local [5] or global optimization methods [6]. The reliability and efficiency
of these methods are successfully being improved by the community, e.g., by
achieving immunity to poor initial designs [7]. Nevertheless, despite reducing
the problem’s complexity, these techniques struggle to reveal new and non-
intuitive designs. Thus, we aim to freely optimize the spatial material
presence in a region representing the design domain to eliminate designer
bias. The task of optimizing material distribution is called topology
optimization [8] which exists in a variety of formulations.

Topology optimization based on genetic algorithms [9] and exact reanalysis
[10, 11], attempts to confront the underlying binary nature of the problem
[12] to find an approximate solution in a large solution space with many local
minima. However, these methods are often limited to problems with relatively
few design variables due to the necessity of solving a state equation one more
time for each design variable for each iteration, resulting in huge amounts of
function evaluations [13].

In contrast to the aforementioned combinatorial methods, density-based
topology optimization [8] applies a continuous relaxation to the binary
problem, i.e., it replaces binary design variables with continuous ones, and,
thus, can exploit adjoint sensitivity analysis [14] to compute sensitivities
independently on the number of design variables. Consequently, gradient
information for objectives is evaluated without significant computational
overhead since only one extra adjoint equation needs to be solved per
iteration. Density-based topology optimization is a deterministic local
algorithm that can converge within a few hundred function evaluations for even
a billion design variables [15].

The density formulation for topology optimization has been intensely
researched during recent decades and has been successfully utilized in various
scientific and engineering branches. Apart from the original mechanical
problems [16], [8], it has also been utilized in electromagnetic applications
with various numerical techniques, such as the optimization of the
distribution of dielectric material, based on finite element modeling (FEM),
e.g., to confine light in photonic crystals [17], and to enhance the
transmission properties of nano-photonic waveguides [18], or, based on finite-
difference time-domain modeling (FDTD), e.g., to improve the bandwidth of
dielectric resonator antennas [19]. A first approach distributing conductive
material was introduced in [20] and extensively studied in FEM to optimize
$s$-parameters of microwave waveguide filters [21] to maximize the radiation
efficiency of electrically small antennas [20], or, based on FDTD, to minimize
return loss of monopole antennas [22].

Density-based topology optimization has also been combined with method of
moments (MoM) for the maximization of the total efficiency of an antenna [23].
In that work, the authors did not regularize the small structural features by
density filtering [24]; thus, their designs are mesh-dependent. Furthermore,
the fractional bandwidth was not considered as an objective. In [25], the
bandwidth is optimized through the associated (inversely proportional, [26])
Q-factor by topology optimization based on genetic algorithms [9] producing
highly irregular designs.

This work aims to develop an automated procedure for optimizing parameters of
electrically small antennas [2] within the MoM paradigm [27]. In particular,
the Q-factor [26] is minimized in this paper to demonstrate the efficacy of
the procedure.

The projection approach [28] is employed in this paper for gradient-based
topology optimization, with the density filtering approach acting as an
intermediate step to converge to a similar design with mesh refinement [24].
The projection scheme is utilized to remedy the introduction of intermediate
densities, known as gray scales, which are difficult to interpret in terms of
realistic material properties and to ensure convergence to near-binary
results. If a conversion from a near-binary to a pure-binary design,
consisting only of void and material, is performed, only a slight drop in
performance is expected.

The performance of the proposed technique is assessed based on distance to
fundamental bounds [29]. It is demonstrated on two examples of design domains
where human-generated designs exist and possess low values of Q-factor, i.e.,
a rectangular region and a meanderline [30] or a spherical region and a
spherical helix [31], [32]. The efficiency of these empirical designs can be
confirmed by the proposed technique delivering similar topologies or refuted
by providing novel better-performing structures.

The paper is structured as follows. The basics of density-based topology
optimization are introduced into MoM models in Section II, and the design
parametrization is proposed. A suitable interpolation function associating
intermediate values of design variables to material parameters is also
discussed. Filtering techniques are developed and utilized to penalize
unwanted behavior in the designs. A flowchart of the computations of the
developed topology optimization framework is summarized in Section IV. A code
for a density-based topology optimization minimizing Q-factor, with annotated
comments, is downloadable at [33]. Section V proposes a method that minimizes
the Q-factor on the two design domains. The method’s technicalities are
discussed and its time performance is compared to topology optimization by
genetic algorithms [9] and to the method based on exact reanalysis [10]. The
quality of existing empirical designs on a spherical shell is demonstrated by
the proposed technique and briefly discussed. The paper is concluded in
Section VI.

##  II Topology Optimization in Method-of-Moments Models

Highly conductive surfaces are considered for optimization throughout this
paper. The design region is delimited by a bounding box in which the optimized
structure is to be located, see the red dashed line in Fig. 1A. The
electromagnetic field can be unambiguously represented by an equivalent
current density located on the design surface and by the electric field
integral equations (EFIE) [34], see Appendix A for further details. The
physical problem is made mathematically tractable by discretizing the design
domain where unknown state variables will reside. Consequently, current
density is expanded into a weighted sum of basis functions and Galerkin’s
procedure [27] is applied to solve the system equation. Without loss of
generality, let us consider triangular discretization and expansion of current
density into a set of $N$ Rao-Wilton-Glisson basis functions [35]. This leads
to a system of linear equations for unknown current expansion coefficients
$\mathbf{I}\in\mathbb{C}^{N\times 1}$ reading

$\mathbf{ZI}=\left(\mathbf{Z}_{0}+\mathbf{R}_{\rho}\right)\mathbf{I}=\mathbf{V},$
(1)

where $\mathbf{Z}_{0}\in\mathbb{C}^{N\times N}$ and
$\mathbf{R}_{\rho}\in\mathbb{R}^{N\times N}$ are the vacuum and material parts
of the impedance matrix, respectively. Vector
$\mathbf{V}\in\mathbb{C}^{N\times 1}$ represents excitation coefficients,
e.g., a plane-wave or a delta gap feeding. Throughout this paper, the
excitation vector $\mathbf{V}$ is given a priori and is, therefore, fixed
during the optimization. The optimization of volumetric dielectric bodies can
be equally defined with minimum modifications employing a volumetric electric
field integral equation, however, such an approach is out of the scope of this
paper.

For a given excitation $\mathbf{V}$, the system (1) defines a physically
admissible response (current density)

$\mathbf{I}=(\mathbf{Z}_{0}+\mathbf{R}_{\rho})^{-1}\mathbf{V},$ (2)

and, therefore, formula (1) is a fundamental optimization constraint.

Figure 1: (A) Bounding box for the rectangular design domain of $L\times L/2$
dimensions, (B) design domain discretized into $T$ triangles with a particular
binary layout of material represented by material function $\chi_{t}$, (C)
continuous distribution of material represented by variable $\rho_{t}$.
Variable $\rho_{t}$ represents a density in each element of the grid.

###  II-A Design Parametrization

In this paper, the goal of topology optimization is to distribute a highly
conducting material, a perfect electric conductor (PEC), over a fixed design
domain to meet the requirements of antenna performance. To facilitate
optimization, the design domain is the same as the modeling domain, and each
triangular element, referred to as a “pixel”, is parametrized with surface
resistivity $R_{\mathrm{s}}$. Hence, the topology optimization task is reduced
to a classical material distribution task where we aim to determine a material
characteristic function

$\chi_{t}(\boldsymbol{r})=\begin{cases}0\quad\iff\quad
R_{\mathrm{s}}\approx\infty\;\Omega~{}\mathrm{(vacuum)},\\\ 1\quad\iff\quad
R_{\mathrm{s}}\approx 0\;\Omega~{}\mathrm{(PEC)},\end{cases}$ (3)

where $\chi_{t}$ controls the layout of two materials: air and a PEC, see Fig.
1B.

Density-based topology optimization [8] commonly replaces the discrete
material function $\chi_{t}$ with the continuous design variable as

$0\leq\rho_{t}\leq 1,$ (4)

so that the design variable is allowed to take intermediate values between the
air and the PEC representing a “density” of material, see Fig. 1C. Subscript t
represents design variable $\rho_{t}$ as a piece-wise constant over the
triangular domain given by the $t$-th triangle.

Consequently, the interpolation function [36] is used to associate surface
resistivity $R_{\mathrm{s}}$ with intermediate values of $\rho_{t}$ and the
physics is projected to material matrix $\mathbf{R}_{\rho}$ as

$\mathbf{R}_{\rho}=\displaystyle\sum_{t=1}^{T}R_{\mathrm{s}}(\rho_{t})\mathbf{\Psi}_{\rho}^{t},$
(5)

where each pixel attains real constant surface resistivity $R_{\mathrm{s}}$
parametrized by design variable $\rho_{t}$ through an interpolation function,
and $\mathbf{\Psi}_{\rho}^{t}$ is a material element matrix, see Appendix A
for its form.

Generating a suitable type of interpolation function is, generally, the most
crucial step within topology optimization algorithms [36]. The function should
vary monotonically, and a slight change in $\rho_{t}$ should produce a non-
negligible response in function value $R_{\mathrm{s}}$. Surface resistivity
acts as a damping parameter for electromagnetic field values, i.e., its
magnitude has a strong effect, and thus, the range of resistivities in the
design plays a crucial role, both for the physics evaluation and optimization
process. Consequently, we cannot adapt the simple linear interpolation
function, which has been successfully employed in dielectric designs [37],
[18], [19], mainly because of the huge contrast between the surface
resistivity of vacuum and PEC. Therefore, a more suitable interpolation scheme
for our design approach is developed in Section II-B.

###  II-B Interpolation Scheme

Considering the aforementioned requirements on the interpolation function, it
is appropriate to express surface resistivity according to [20] and [37] as

$R_{\mathrm{s}}(\rho_{t})=R_{\mathrm{vac}}\left(\dfrac{R_{\mathrm{met}}}{R_{\mathrm{vac}}}\right)^{f(\rho_{t})},$
(6)

where $f(\rho_{t})$ is a penalization function based on a RAMP-like111Rational
Approximation of Material Properties. model [8]

$f(\rho_{t})=\frac{\rho_{t}}{1+p(1-\rho_{t})},$ (7)

where $p$ is a penalty parameter used to regulate the behavior of the
interpolation and penalize intermediate values. In this work, $p=1$, and we
utilize other penalization techniques shown later in Section II-D. Surface
resistivity (6) attains range $[R_{\mathrm{vac}},R_{\mathrm{met}}]$, whose
boundaries represent the resistivity of the vacuum and PEC.

The upper bound $R_{\mathrm{met}}=R_{\mathrm{s}}(\rho_{t}=1)$ and the lower
bound $R_{\mathrm{vac}}=R_{\mathrm{s}}(\rho_{t}=0)$ on attainable
resistivities should be set to simulate the actual physical response of the
material distribution correctly but, at the same time, have the smallest
possible contrast. They are set based on the numerical investigations of the
authors, see Appendix B, as

$R_{\mathrm{vac}}=10^{5}\;\Omega,\quad R_{\mathrm{met}}=1\;\Omega,$ (8)

which, for this paper (Q-factor minimization), provides physically acceptable
results and enough freedom for the optimization step since the Q-factor is not
influenced by loss power, see Appendix B for more details. Therefore,
resistivity bounds (8) represent the approximative physical behavior of the
vacuum and PEC, but a subsequent analysis of the final design using the true
resistivities, i.e., $R_{\mathrm{vac}}\approx\infty\;\Omega$ and
$R_{\mathrm{met}}\approx 0\;\Omega$, results in a negligible drop in Q-factor
value. Nevertheless, current density is influenced, yet it is important to
emphasize that these values are not universal and cannot be used for arbitrary
physical metrics. Numerical investigations should always be conducted to
determine the properties for a particular optimization goal.

The interpolation function (6) is shown in Fig. 2 for three values of
$R_{\mathrm{vac}}$. The exponential dependency on $\rho_{t}$ results in a fast
decay from $\rho_{t}=1$, which is demonstrated on a current distribution on a
center-fed dipole antenna in which the end parts are covered by variable
resistivity ($\rho_{t}$ is varied from 0 to 1) while the central part is
covered with resistivity $R_{\mathrm{met}}=R_{\mathrm{s}}(\rho_{t}=1)$. We
observe standard cosinusoidal and triangular current distribution for the
half-wave and a short dipole, respectively.

Figure 2: Interpolation scheme of surface resistivity
$R_{\mathrm{s}}(\rho_{t})$ as a function of density $\rho_{t}$ for three
different values of vacuum resistivity $R_{\mathrm{vac}}$. The influence of
the interpolation on a current density is explored on a half-wave-long strip
dipole antenna with the ends being designable domain $\varOmega_{\mathrm{d}}$,
in which $\rho_{t}$ is varied from 0 to 1. The Red solid lines show a current
distribution along the dipole excited by a delta gap placed in the middle,
e.g., a standard $\lambda/2$-long dipole current is observed for $\rho_{t}=1$.
Figure 3: Projection filter $\mathcal{H}$ with four different values of
sharpness $\beta$ and projection level $\eta=0.5$. Figure 4:  Demonstration of
the filtering process. The initial test design (solid box) represented by
$\rho_{t}$ contains specific patterns, e.g., sheets, and holes with various
widths or checkerboard patterns, comparable in size with density filter radius
$R_{\mathrm{min}}$. Density filtering of these patterns results in
distribution $\widetilde{\rho_{t}}$ containing gray areas (dashed boxes),
possibly removing undesirable details by smoothing the design distribution,
thus blurring the solution space which might result in an escape from the
local minimum. The subsequent application of the projection filter removes
unwanted gray areas resulting in a design of near-binary nature (dotted
boxes).

###  II-C Minimization of Q-Factor

For single resonance systems, Q-factor provides a measure of the frequency
selectivity of an antenna [38], i.e., an estimate of a fractional bandwidth
through its inverse proportionality $B\sim Q^{-1}$ [26], see Appendix C for
its definition. For electrically small antennas the lowest Q-factor is
observed at self-resonance [29], [39]. Nevertheless, previous design
experience is required to construct such a structure, e.g., a spherical helix
antenna [32], or a meanderline antenna [30].

The density-based topology optimization is carried out in this paper to find
an optimized layout of conductive material to minimize the antenna Q-factor
within MoM models. The optimization formulation is given by

$\displaystyle\underset{\rho_{t}}{\mathrm{minimize}}$ $\displaystyle\quad
Q(\mathbf{I})=\mathrm{max}\\{Q_{\mathrm{e}}(\mathbf{I}),Q_{\mathrm{m}}(\mathbf{I})\\}$
(9) $\displaystyle\mathrm{subject~{}to}$
$\displaystyle\quad\left[\mathbf{Z}_{0}+\mathbf{R}_{\rho}(\rho_{t})\right]\mathbf{I}=\mathbf{V},$
$\displaystyle\quad\dfrac{1}{A_{0}}\sum_{t=1}^{T}\rho_{t}A_{t}\leq S_{f},$
$\displaystyle\quad 0\leq\rho_{t}\leq 1,\quad t=1,\dotsc,T,$

where subscripts e/m denote the electric and magnetic Q-factor parts,
respectively, $\rho_{t}$ is the density design variable associated with each
triangular pixel, $A_{t}$ is the area of the $t$-th triangle, $A_{0}$ is the
design domain area, and $S_{f}$ is the maximum fraction of design domain area
which can be occupied with material. Task (9) is defined as a nested
formulation, i.e., the MoM equation is included for the sake of completeness
but its solution is directly incorporated into the evaluation of the objective
and constraints [8]. Although electromagnetic designs do not benefit from
spanning the whole design domain with material, the second constraint, which
restricts the amount of the available material, is applied to speed up the
convergence and to limit the existence of non-connected regions of the
conductor.

###  II-D Filtering Techniques

Solving topology optimization task (9) as it stands often results in mesh-
dependent designs. Whenever a mesh is refined, the topology optimization may
add details comparable to the mesh element size and often converge to a
significantly different topology. Moreover, checkerboard-like patterns might
occur in this resulting design [40]. The continuous relaxation also inevitably
leads to the appearance of residual gray pixels, i.e., intermediate densities,
in the actual design. This results in final designs being practically
unrealizable.

A projection operator [28] is used to project continuous design variables into
a binary space as

${\bar{\widetilde{\rho_{t}}}}=\mathcal{H}({{\widetilde{\rho_{t}}}}),$ (10)

where $\bar{\widetilde{\rho_{t}}}$ is the projected filtered density and
$\mathcal{H}$ is the Heaviside function. The Heaviside function does not have
a finite derivative and is thus replaced by a continuous function [28] in the
form

$\mathcal{H}({\widetilde{\rho_{t}}})\approx\frac{\tanh(\beta\eta)+\tanh(\beta({\widetilde{\rho_{t}}}-\eta))}{\tanh(\beta\eta)+\tanh(\beta(1-\eta))},$
(11)

where parameters $\beta$ and $\eta$ are threshold sharpness and threshold
level, respectively, see Fig. 3.

The original design variables $\rho_{t}$ are smoothed to regularize the
solution space of the optimization problem to provide an existence of
solutions. Several regularization schemes are reported in the literature,
e.g., perimeter control [41], sensitivity filtering [42], and density
filtering [24], [43]. In this paper, density filtering operation
$\mathcal{D}(\cdot)$ is used as an intermediate step written in the form of a
convolution [43] as

$\widetilde{\rho_{t}}=\mathcal{D}(\rho_{t})=\sum_{j\in B_{t}}W_{tj}\rho_{j},$
(12)

where $\widetilde{\rho_{t}}$ is the filtered variable, $\rho_{j}$ is the
original design variable and $W_{tj}$ is a convolution kernel defined as

$W_{tj}=\frac{(R_{\mathrm{min}}-\|\boldsymbol{r}_{t}-\boldsymbol{r}_{j}\|)}{\sum\limits_{{k}\in
B_{{t}}}(R_{\mathrm{min}}-\|\boldsymbol{r}_{t}-\boldsymbol{r}_{k}\|)},$ (13)

with $\boldsymbol{r}_{m}$ being the center of the $m$-th triangle. Support
domain $B_{t}$ for the filter is a mesh-independent neighborhood of $t$-th
triangle specified by a sphere with radius $R_{\mathrm{min}}$ as

$B_{t}\in\left\\{j\;:\;\|\boldsymbol{r}_{t}-\boldsymbol{r}_{j}\|\leq
R_{\mathrm{min}}\right\\}.$ (14)

The density filtering is introduced to impose a length-scale, i.e., a
characteristic length used to control the size of features in the optimized
design. It can be interpreted as the minimum feature size the optimizer can
generate. The density filter is utilized to regularize the solution space and
accelerate the optimization process but introduces intermediate densities ,
motivating the subsequent projection towards binary design by (11).

The projection operator (11) is utilized in a continuation scheme to converge
to a near-binary design in order to moderate adverse effects. Sharpness
$\beta$ is gradually increased during the optimization to minimize the effect
of gray elements and to converge to a near-binary result. The pure-binary
material distribution is obtained by the hard thresholding so that
$\bar{\widetilde{\rho_{t}}}\geq 0.5\to 1$ and
$\bar{\widetilde{\rho_{t}}}<0.5\to 0$, i.e., applying (11) with sharpness
$\beta\to\infty$.

The influence of filters is demonstrated in Fig. 4, where a test design with
patterns of varying complexity is filtered by a density filter with two
different radii $R_{\mathrm{min}}$ comparable with the included details.
Despite the introduction of gray areas by the density filter, unwanted
details, e.g., the checkerboard pattern, are removed. The density filtering
technique provides the existence of the solution and convergence with the mesh
refinement. The subsequent utilization of the projection filter (11) limits
the gray areas and ensures convergence to near-binary results, but the length-
scale imposed by the density filter is lost. To ensure length-scale control as
well, a robust formulation of topology optimization [28] may be exploited, but
such a formulation is out of the scope of this work.

##  III Q-factor Minimization via Topology Optimization

The optimization task (9) is regularized by the filtering techniques, i.e.,
density, and projection filters, and requires continuous derivatives to be
solved by the gradient-based optimizer. In order to solve the originally non-
differentiable “min-max” problem (9), the problem is converted to the so-
called bound formulation [8]

$\displaystyle\underset{\rho_{t}}{\mathrm{minimize}}$ $\displaystyle\quad z$
(15) $\displaystyle\mathrm{subject~{}to}$ $\displaystyle\quad z\geq
Q_{\mathrm{e}}(\mathbf{I}),$ $\displaystyle\quad z\geq
Q_{\mathrm{m}}(\mathbf{I}),$
$\displaystyle\quad\left[\mathbf{Z}_{0}+\mathbf{R}_{\rho}(\bar{\widetilde{\rho_{t}}})\right]\mathbf{I}=\mathbf{V},$
$\displaystyle\quad\dfrac{1}{A_{0}}\sum_{t=1}^{T}\bar{\widetilde{\rho_{t}}}A_{t}\leq
S_{f},$
$\displaystyle\quad\bar{\widetilde{\rho_{t}}}=\mathcal{H}\left(\mathcal{D}\left(\rho_{t}\right)\right),$
$\displaystyle\quad 0\leq\rho_{t}\leq 1,\quad t=1,\dotsc,T,$

where an extra decision variable $z$ is introduced. We caution, however, that
the optimization task (LABEL:eq:MMA_task) is continuous but non-convex, thus
global optimality is not guaranteed. Therefore, the quality of the solution
can be solely judged based on the comparison to fundamental bounds [29] or
comparison to existing designs.

The topology optimization task (LABEL:eq:MMA_task) is solved for local
optimality using the Method of Moving Asymptotes (MMA) [44], which is a
gradient-based optimizer requiring knowledge of the sensitivities of
objectives and constraints with respect to the design variable. We here use
that all components of objectives and constraints are continuously
differentiable. The sensitivities of $Q_{\mathrm{e}}$ and $Q_{\mathrm{m}}$,
with respect to physical design variable $\bar{\widetilde{\rho_{t}}}$, are
evaluated by an adjoint sensitivity analysis [14] as

$\frac{\,\mathrm{d}Q_{\mathrm{e/m}}}{\,\mathrm{d}\bar{\widetilde{\rho_{t}}}}=2\mathrm{Re}\left\\{\boldsymbol{\lambda}^{\mathrm{T}}_{\mathrm{e/m}}\frac{\partial\mathbf{R}_{\rho}}{\partial\bar{\widetilde{\rho_{t}}}}\mathbf{I}\right\\},$
(16)

where adjoint vectors $\boldsymbol{\lambda}_{\mathrm{e/m}}$ are determined
from the corresponding adjoint equations

$\mathbf{Z}^{\mathrm{T}}\boldsymbol{\lambda}_{\mathrm{e/m}}=-\left(\frac{\partial
Q_{\mathrm{e/m}}}{\partial\mathbf{I}}\right)^{\mathrm{T}}.$ (17)

The derivative of $\mathbf{R}_{\rho}$ with respect to
$\bar{\widetilde{\rho_{t}}}$ is

$\frac{\partial\mathbf{R}_{\rho}}{\partial\bar{\widetilde{\rho_{t}}}}=\frac{\partial
R_{\mathrm{s}}}{\partial\bar{\widetilde{\rho_{t}}}}\mathbf{\Psi}_{\rho}^{t},$
(18)

where $\partial R_{\mathrm{s}}/\partial\bar{\widetilde{\rho_{t}}}$ is the
derivative of the interpolation function (6). Thus, to compute the sensitivity
of the Q-factor, only the derivatives $\partial
Q_{\mathrm{e}}/\partial\mathbf{I}$ and $\partial
Q_{\mathrm{m}}/\partial\mathbf{I}$, see Appendix C, need to be evaluated. The
system matrix is reused to solve adjoint equations (17) with the corresponding
right-hand side.

Sensitivities (16) are evaluated with respect to the filtered and projected,
i.e., physical, density field $\bar{\widetilde{\rho_{t}}}$, hence, the chain
rule is utilized to determine sensitivities with respect to a change in the
design field $\rho_{t}$ as

$\frac{\,\mathrm{d}Q_{\mathrm{e/m}}}{\,\mathrm{d}\rho_{t}}=\sum_{i\in
B_{t}}\frac{\,\mathrm{d}Q_{\mathrm{e/m}}}{\,\mathrm{d}\bar{\widetilde{\rho_{i}}}}\frac{\partial\bar{\widetilde{\rho_{i}}}}{\partial\widetilde{\rho_{i}}}\frac{\partial\widetilde{\rho_{i}}}{\partial\rho_{t}},$
(19)

where the derivative of the projected density $\bar{\widetilde{\rho_{t}}}$,
with respect to a change in the filtered density $\widetilde{\rho_{t}}$, is
determined as

$\frac{\partial\bar{\widetilde{\rho_{i}}}}{\partial\widetilde{\rho_{i}}}=\frac{\beta(1-\tanh(\beta(\widetilde{\rho_{i}}-\eta))^{2}}{\tanh(\beta\eta)+\tanh(\beta(1-\eta))},$
(20)

and where the sensitivity of the filtered field $\widetilde{\rho_{t}}$, with
respect to design variable $\rho_{t}$, is found as

$\frac{\partial\widetilde{\rho_{i}}}{\partial\rho_{t}}=\frac{(R_{\mathrm{min}}-\|\boldsymbol{r}_{t}-\boldsymbol{r}_{i}\|)}{\sum\limits_{{k}\in
B_{t}}(R_{\mathrm{min}}-\|\boldsymbol{r}_{t}-\boldsymbol{r}_{{k}}\|)}.$ (21)

Figure 5: Flowchart of the proposed topology optimization procedure.

##  IV Flowchart of the Topology Optimization procedure

To implement the optimization procedure, users can refer to the code
implementation of density topology optimization which is available as
supplementary material [33] and includes annotated comments for each part of
the process. The flowchart in Fig. 5 provides an overview of the optimization
stages and is further detailed in this section. During the preprocessing, MoM
operators are calculated on a fixed grid, and the user input parameters, such
as area fraction ($S_{f}$), density filter radius ($R_{\mathrm{min}}$),
initial design variable ($\rho_{t}$), maximum iterations ($I$), and maximum
design variable change ($\Delta\rho_{t}^{\mathrm{max}}$, defaulted to 0.01)
are set. The projection filter (11) starts at level $\eta=0.5$ and sharpness
$\beta=1$, which is later updated in the optimization loop until it reaches a
maximum of $\beta_{\mathrm{max}}=32$.

The iterative optimization process follows: filtering, solving the MoM
equation, evaluating objectives, computing sensitivities, and updating the
design variable ($\rho_{t}$) by MMA with move limits [44] set to 0.25. If the
maximum number of iterations is not reached and the maximum change in
$\rho_{t}$ in two consecutive iterations $i$ is greater than
$\Delta\rho_{t}^{\mathrm{max}}$, i.e., if conditions

$\underset{{t}}{\mathrm{max}}\big{\\{}|\rho_{t}^{i}-\rho_{t}^{i-1}|\big{\\}}\geq\Delta\rho_{t}^{\mathrm{max}}\quad\mathrm{and}\quad
i<I$ (22)

are met, the optimization continues, and a decision about the update of
sharpness $\beta$ of the projection filter is made. The sharpness is doubled
if the change in $\rho_{t}$ in two consecutive iterations $i$ is less than
$\Delta\rho_{t}^{\mathrm{max}}$ or every 100 iterations until the sharpness
reaches maximum value, i.e., if logic relation

$\displaystyle\big{(}\underset{{t}}{\mathrm{max}}\big{\\{}|\rho_{t}^{i}-\rho_{t}^{i-1}|\big{\\}}\leq\Delta\rho_{t}^{\mathrm{max}}$
$\displaystyle\mathrm{or}\quad\mathrm{mod}(i,100)=1\big{)}$ (23)
$\displaystyle\mathrm{and}\quad\beta\leq\beta_{\mathrm{max}}$

are fulfilled, and where $\mathrm{mod}()$ is the modulo operation. After
meeting the stopping criterion, the design and physical quantities are post-
processed and visualized. The continuation scheme for the projection filter
produces a near-binary optimized design which is finally converted to a pure-
binary design through hard thresholding. The Q-factor of the thresholded
design is denoted as $\raisebox{1.5pt}{\hbox to0.0pt{$\ulcorner$\hss}\kern
2.0pt}Q$ and is evaluated with physically accurate values of the resistivities
for the vacuum ($R_{\mathrm{vac}}\approx\infty\;\Omega$) and PEC
($R_{\mathrm{met}}\approx 0\;\Omega$).

Optimization parameters significantly impact the quality of the solutions, and
multiple optimization runs with different initial guesses or parameters may be
necessary to reach the best possible optimized design as the problem is non-
convex and can become trapped in local minima. Care must also be taken in
selecting $R_{\mathrm{min}}$ and the continuation scheme’s speed, as a large
filter radius can prevent good minima from being reached, and a rapid increase
in $\beta$ can destabilize convergence. Testing different parameters on a
coarse mesh before running the final optimization with a fine grid is highly
recommended.

Figure 6: (Top) Convergence history of Q-factors
$Q=\mathrm{max}\\{Q_{\mathrm{e}},Q_{\mathrm{m}}\\}$ normalized to the
fundamental bound $Q_{\mathrm{lb}}^{\mathrm{TM}}$ [45], where dashed vertical
lines represent iterations in which the sharpness $\beta$ of projection filter
$\mathcal{H}$ was doubled. (Bottom) Snapshots of density field distribution
from a few selected iterations are included below in the convergence graph.
The optimizer reaches a self-resonant design in iteration $i=29$. The
structure is only slightly refined during the rest of the iterations, mainly
by removing gray elements by the continuation scheme of the projection filter.
The black line bounds the region where $\bar{\widetilde{\rho_{t}}}\geq 0.01$.
The density filter size characterized by radius $R_{\mathrm{min}}$ is also
depicted in scale with the consequent utilization of the projection filter
with sharpness $\beta$, effectively reducing its size and losing the imposed
length-scale.

##  V Numerical Investigations

The proposed topology optimization technique is implemented in MATLAB [46]
according to Section IV, and its properties and performance are discussed. The
goal is to minimize the Q-factor (LABEL:eq:MMA_task) by distributing a
conductive material over a fixed triangular discretized domain. The required
MoM matrices are computed in AToM [47]. The examples are evaluated on a
computer with Intel Xeon Gold 6244 CPU (16 physical cores, 3.6 GHz) with 384
GB RAM. The excitation of the antennas is performed via delta gap feeders with
their immediate neighborhood kept fixed222The excitation must be connected at
all times, and the action of the density filter may cause unwanted
disconnections or shortages of the feed..

###  V-A Rectangular Design Region

The first example considers a rectangular design region of aspect ratio $3:5$
and electrical size $ka=0.8$, where $k$ is the wavenumber and $a$ is the
radius of the smallest circumscribing sphere. A delta gap feeder provides
excitation close to the left boundary, see the design at the bottom of Fig. 6.
The objectives, i.e., Q-factors $Q_{\mathrm{e/m}}$, are normalized to the
bound $Q_{\mathrm{lb}}^{\mathrm{TM}}$ [45] which is valid for antennas
radiating TM modes only. The design domain is discretized into 3 840 triangles
(5 696 basis functions). The Q-factor minimization (LABEL:eq:MMA_task) is
performed by the density-based topology optimization with a fixed radius of
the density filter set to $R_{\mathrm{min}}=0.15a$, gradually increased
sharpness of the projection filter according to Section IV, and the maximal
number of iterations $I=600$. The optimization is subjected to area fraction
constraint $S_{f}=0.35$, and all design variables $\rho_{t}$ are initially set
to the prescribed area fraction.

The optimization ran for 135 minutes, and its progress is reported at the top
of Fig. 6. It can be observed that the proposed technique quickly locates a
self-resonant structure that separates the electric charge well, i.e., an
expected solution of the minimal Q-factor [48]. The subsequent convergence
rate is mostly flat, with occasional jumps resulting from the projection
filter update. The design is only slightly refined by further removing the
gray elements. The optimized design, see $i=450$ at the bottom of Fig. 6,
exhibits a Q-factor very close to the fundamental bound
$Q/Q_{\mathrm{lb}}^{\mathrm{TM}}=1.06$. Below the snapshots in Fig. 6, you
will find information about the impact of the projection filter on the density
filter size while updating sharpness $\beta$. While the projection scheme
yields a near-binary design, it also disrupts the length-scale imposed by the
underlying density filter. If no special precautions are taken, the length-
scale near the outer boundary may be violated [49]. The exact length-scale in
this paper is not crucial as long as the final design avoids unwanted
features, e.g., checkerboard-like patterns [40].The primary purpose of
filtering is to regularize the solution space and accelerate the optimization
process.

Figure 7: Performance of topology optimization algorithms for Q-factor
minimization for two grids. Settings are the same as in Fig. 6. The objective
is normalized to the fundamental bound $Q_{\mathrm{lb}}^{\mathrm{TM}}$. Three
implementations are shown: density topology optimization (black curves), a
genetic algorithm (green curves) and a memetic combination of global and local
step (red curves). Black crosses mark the computational time needed to reach
$Q=2Q_{\mathrm{lb}}^{\mathrm{TM}}$ and show the rough estimate on the
scalability of the method with the number of optimization variables.

Despite the proposed method being an inherently local procedure, it is
resilient towards initial design choices due to the nature of the Q-factor
[50] and the usage of spatial filtering, see Appendix D for more details.
While optimizing the Q-factor, the method delivers locally extremal designs
that are in the performance close to the fundamental bound, thus approaching a
global extremum.

The computational time necessary to obtain a design from the proposed
optimization is further judged for two discretization grids, see Fig. 7. It is
observed that the density topology optimization provides convergence to local
minimum and has a steep convergence rate which is typical for gradient-based
algorithms. It provides gray designs with low values of the Q-factor in a
short computational time. Topology optimization using density-based approaches
exhibits good scalability with increasing design variables, as noted in [15].
For a rough estimate, please refer to the gray region in Fig. 7. The
relatively low computational cost allows rerunning the optimization several
times for various preferences of user-defined parameters, e.g., interpolation
function, area fraction, etc., to determine which settings produce the best
result.

The computation time is also compared with other topology optimization
techniques, namely, binary topology optimization based on a genetic algorithm
and topology optimization based on memetic scheme [10] combining exact
reanalysis and a genetic algorithm, see Fig. 7. Both methods utilize 144
agents. For the sake of simplicity, the computation time is not normalized to
the number of physical cores used for the evaluation. The convergence rate of
the memetic scheme is superior to the sole use of the genetic algorithm, but
it is much slower than the density-based approach. An estimate on scalability
is also included, see the red and green regions in Fig. 7, and is inferior to
the density topology optimization.

Figure 8: (Top) Frequency sweep of the normalized Q-factor of the design
optimized at single electrical size $ka=0.8$. The purple curve corresponds to
density-based topology optimization with gradually increased sharpness of the
projection filter. The orange curve corresponds to the thresholded design. The
thresholding may result in a slight frequency shift and performance drop.
(Bottom left) Optimized and (bottom right) thresholded designs are also
included with the dipole-like current density excited by a delta gap feeder
depicted by the blue line. The colormap shows an absolute value of the current
density. The black line bounds the region where
$\bar{\widetilde{\rho_{t}}}\geq 0.01$. Figure 9: The proposed optimization
technique is utilized to minimize the Q-factor at the various electric sizes
and fixed discretization ($T=3840$). The first row shows the thresholded
design and corresponding Q-factor. The grid cannot accommodate a self-resonant
structure for the small electric sizes, and much finer granularity would be
required. In the second row, the optimization is repeated for a particular
electric size $ka=0.4$ while the discretization grid is changed to an
unstructured one and further refined. The density filter with a fixed radius
provides the solution while the mesh is refined.

Despite the use of a continuation scheme, the optimized design obtained by the
density topology optimization still contains residual gray elements, see the
design in the bottom left corner of Fig. 8. The pure-binary design obtained by
the hard thresholding is included in the bottom right corner of Fig. 8. The
consequences of the hard thresholding, e.g., performance drop, and resonance
shift [51], [52], are moderated by the continuation scheme but not entirely
suppressed. The binary design still exhibits low Q-factor
$\raisebox{1.5pt}{\hbox to0.0pt{$\ulcorner$\hss}\kern
2.0pt}Q/{Q_{\mathrm{lb}}^{\mathrm{TM}}}=1.06$, but the self-resonance is
slightly shifted, see the frequency sweep in Fig. 8.

We repeatedly perform Q-factor minimization on the identical grid ($T=3840$)
for $ka={0.1,0.4,0.8,1.2}$ to demonstrate the evolution of the optimized
designs with electrical size. The designs are hard thresholded and shown in
the first row of Fig. 9. It can be observed that a finer grid is required for
low electric sizes $ka=\\{0.1,0.4\\}$ to accommodate a structure with a self-
resonant Q-factor, e.g., a meander line design [30], to reach the lower bound.
The optimizer realizes a self-resonant structure at $ka=1.2$, but it perhaps
requires a different feed position to get the lower bound.

To demonstrate the effect of the density filter, we repeat the optimization
for a fixed electric size of $ka=0.4$ while refining the discretization grid.
Please refer to the thresholded designs in the second row of Fig. 9. The
density filter with fixed radius $R_{\mathrm{min}}$ regularizes the solution
space and can achieve structural convergence while the mesh is refined, see
similarities in the designs. Nevertheless, the precise length-scale imposed by
the density filter is lost while utilizing the continuation scheme of the
projection filter [28].

By visual inspection, the mesh with the finest granularity can accommodate a
self-resonant meander line antenna; however, the density filter size must be
further reduced, which may destabilize the convergence.

Figure 10: Comparison of empirical and optimized designs, in terms of Q-factor
values, are presented. The red line depicts self-resonant 2-arm spherical
helix antennas, and the green line shows self-resonant 2-arm loxodrome
antennas designed for a given electrical size $ka$. Purple circles represent a
single topology optimization run for the particular settings of
$R_{\mathrm{min}}$ and $S_{f}$, i.e., ten runs are performed for each $ka$. A
purple line connects the best optimized designs, and the orange crosses depict
the best thresholded structures. Notice that the discretization grid cannot
accommodate a 2-arm geometry for low values of $ka$. For that reason, the
topology optimization realizes a self-resonant 1-arm structure.

###  V-B Spherical Design Domain

In the second example, the density topology optimization is utilized to solve
the optimization problem (LABEL:eq:MMA_task) for a spherical shell discretized
into 2 400 triangular elements (3 600 basis function). The structure is
excited by two identical delta gap sources placed symmetrically on the
opposite sides of the equator to cover a shell uniformly with the material.
The performance of the optimized designs is assessed by the fundamental bound
[53], which for small electrical sizes reaches $(ka)^{3}Q_{\mathrm{lb}}\approx
1$. The optimization is performed with the maximal number of iterations
$I=600$ for several electrical sizes in an interval $ka\in[0.12,0.28]$,
several radii of the density filter
$R_{\mathrm{min}}=\\{0.16,0.18,0.20,0.22\\}a$, and two area fractions
$S_{f}=\\{0.4,0.5\\}$. As previously mentioned, it is often mandatory to
restart the optimization with various initial parameters due to the non-convex
nature of the problem.

Optimized designs obtained by topology optimization (average computation time
50 minutes) for various parameters are shown in Fig. 10 by the purple circle
marks. The purple line shows the performance of the best designs.
Unsurprisingly, the optimizer finds a structure that combines two dominant
spherical modes [53]. Thresholding the optimized designs results in a slight
drop in performance, see orange line in Fig. 10. For $ka<0.2$, the optimizer
realizes self-resonant 1-arm designs and leaves one delta gap feeder isolated.
A considerably finer mesh grid would be required to accommodate 2-arm geometry
reaching the self-resonance. For $ka>0.2$, the discretization and the selected
radius of the density filter form a solution space in which the proposed
method finds self-resonant 2-arm spherical geometry with a low Q-factor value,
i.e., for $ka=0.24$ $(ka)^{3}Q=1.25$.

State-of-the-art spherical designs are helical antennas [32], see Appendix E
for their parametrization. Their Q-factors can further be reduced by folding,
i.e., using a multi-arm spherical helix, which covers the spherical surface
more uniformly333In [32] the motivation for realizing the folded design was to
increase the radiation resistance. Folding, nevertheless, also helps to
approximate dominant spherical modes better.. Compared with the topology-
optimized designs, Q-factors of empirical 2-arm spherical helices for various
electrical sizes $ka$ are depicted in Fig. 10 by the red dots. The optimized
designs outperform the empirical designs since the optimizer operates with
more design variables. It is also worth noting that the optimized design does
not resemble a spherical helix. The slope of a spherical helix is not constant
along its curve, see Appendix E. A curve conforming with the optimized design
is loxodrome (Rhumb line) [54] which always crosses the meridians at the same
angle, see Appendix E for the parametrization. This agrees with the work [31]
where it was analytically shown that an optimal current path for minimal
Q-factor value follows loxodromes, but the authors improperly called such
curves spherical helices. The authors of this manuscript recommend a
differentiation between these curves since their Q-factor performance is not
the same. Loxodrome is only similar to a spherical helix near the equator but
has an asymptotic point at the poles, i.e., the stereographic projection of a
loxodrome to the plane of the equator is a logarithmic spiral. The Q-factor of
the empirical 2-arm self-resonant loxodrome antennas for electrical sizes
$ka\in[0.12,0.28]$ are depicted in Fig. 10 by the green dots. It can be
observed that loxodromes outperform the 2-arm spherical helix design in the
entire $ka$ region. Curiously, they also outperform the optimized designs for
$ka<0.2$. An insufficient mesh grid density used for the optimization causes
this. For $ka>0.2$, the loxodrome has a higher Q-factor value than the
optimized designs.

Figure 11: (Top) Area fraction $S_{f}$, and (middle) convergence history of
Q-factors $Q=\mathrm{max}\\{Q_{\mathrm{e}},Q_{\mathrm{m}}\\}$ normalized to
the lower bound $Q_{\mathrm{lb}}$ [53] during the optimization. The dashed
vertical lines represent iterations in which the sharpness $\beta$ of
projection filter $\mathcal{H}$ was doubled. (Bottom) Snapshots of density
field distribution from several iterations are included. The black line in the
designs bounds the region where $\bar{\widetilde{\rho_{t}}}\geq 0.5$, thus
showing a thresholded design boundary. Figure 12: The optimized design was
obtained by density topology optimization. The structure’s resemblance to the
two spherical curves, i.e., spherical helix and loxodrome, is depicted. It is
shown that the optimized topology follows the loxodrome curve more closely
than the spherical helix. Both curves are very alike near the equator but
distinct at the poles [31].

The optimization was repeatedly performed for electrical size $ka=0.24$ for a
finer discretization of a spherical shell consisting of 11 616 triangles (17
424 basis functions) to obtain a higher level of structural details. The
feeding configuration remains the same. The fixed radius of the density filter
$R_{\mathrm{min}}=0.2a$ was used, and the optimization was subjected to the
area fraction constraint444On the grid with $T=2400$ the optimizer distributed
material using only 48% of the surface, and thus, $S_{f}$ is reduced to make
the area fraction constraint active. $S_{f}=0.48$. The algorithm ran for 52
hours and successfully provided a design occupying 48% of the design domain,
see the first plot in Fig. 11. The evolution of the optimized Q-factor is
captured in the second plot in Fig. 11, where $Q_{\mathrm{m}}$ and
$Q_{\mathrm{e}}$ approach to realize a self-resonant topology. The snapshots
of the spherical designs during the optimization are included at the bottom of
Fig. 11. Within tens of iterations, the optimizer finds a 2-arm spherical
topology whose arms are further extended by turning around the poles. The
optimized structure has an optimized slope to combine dominant TM and TE
spherical modes properly. The optimized design obtained by the density
topology optimization exhibits $(ka)^{3}Q=1.23$ and is shown in Fig. 12 for
two orthogonal views. The spherical helix and loxodrome curves are also
plotted for comparison.

The analytical approach [31] dictates an optimal current combination of
dominant spherical modes, which forms a constant angle with the meridians. The
loxodrome has the same properties and can approximate the optimal current
distribution better than the spherical helix, see the difference of the curves
mentioned above from the optimized topology in Fig. 12. The optimizer more
accurately follows the loxodrome than the spherical helix. The proposed method
can deliver a close-to-optimal structure in terms of Q-factor value in an
automated manner.

##  VI Conclusion

A complete treatment of density-based topology optimization in a method-of-
moments framework was presented and verified on a particular problem of the
Q-factor minimization utilizing conductor-based electrically small antennas.
The proposed technique utilizes the projection operator in a continuation
scheme to obtain a near-binary and standard density filter as an intermediate
step to smooth an optimized design. The adjoint sensitivity analysis is
exploited to attain gradient information to solve the optimization problem to
local optimality.

The efficacy of the presented approach was demonstrated by applying it to
canonical examples of optimizing the Q-factor of rectangular and spherical
design regions. The method excels in the scalability and efficient evaluation
of gradients compared to other formulations of topology optimization. The
realized performance is, in all cases, close to the fundamental bounds and
self-resonance is met depending on the granularity of the discretization grid.
Furthermore, the proposed approach is capable of improving the performance of
existing empirical designs in this case without modifying the nature of the
topology. Alternatively, density topology optimization can inspire antenna
engineers by providing novel structures used as starting points for future
design methodology or by proving the existence of well-performing design. The
framework can handle multi-port and multi-frequency optimization only at the
cost of the linear increase of computational complexity.

The settings of the proposed method are tailored for Q-factor optimization,
and may not work universally for a general electromagnetic metric. The
optimization parameters significantly impact the optimization process and its
success, and multiple runs with different settings are necessary. While the
method has the advantage of fast convergence, the resulting designs are near-
binary only, and a projection to a pure-binary result is required. Such a
projection may result in a performance drop and resonance shift. Thus, the
proposed technique can be conveniently combined into a two-step approach, in
which the second step would tune the pure-binary design into self-resonance.

There are many possibilities for future research that are closely connected to
the proposed framework. For example, optimizing feeding position, amplitude,
and phase can improve performance further and offer greater versatility. Other
figures of merits should be explored as well. In the future, the framework
could also be extended to include multiobjective optimization. Another
potential direction is to define the topology optimization with a robust
formulation [28] so that it produces designs resilient to manufacturing
errors. A final possibility is to extend the implementation of the method of
moving asymptotes to cover the globally convergent version which is more
effective for highly non-convex systems.

##  Appendix A Electric Field Integral Equations within MoM formalism

Assume a time-harmonic steady state at angular frequency $\omega$ with
convention $\partial/\partial t\to\mathrm{j}\omega$, with $\mathrm{j}$ being
the imaginary unit. Without loss of generality, surface formulation of EFIE
[34] for highly conducting bodies with space-dependent surface resistance
$R_{s}(\boldsymbol{r})$ is employed as

$-\boldsymbol{n}\times\boldsymbol{n}\times\left(\boldsymbol{E}^{\mathrm{i}}(\boldsymbol{r})+\boldsymbol{E}^{\mathrm{s}}(\boldsymbol{r})\right)=R_{s}(\boldsymbol{r})\boldsymbol{J}(\boldsymbol{r}),$
(24)

where $\boldsymbol{n}$ is the surface normal unit vector,
$\boldsymbol{E}^{\mathrm{i}}(\boldsymbol{r})$ is an incident field and
$\boldsymbol{E}^{\mathrm{s}}(\boldsymbol{r})$ is scattered field defined via

$\boldsymbol{E}^{\mathrm{s}}(\boldsymbol{r})=-\mathrm{j}kZ_{0}\int\limits_{\varOmega^{\prime}}\overline{\mathbf{G}}\left(\boldsymbol{r},\boldsymbol{r^{\prime}}\right)\cdot\boldsymbol{J}\left(\boldsymbol{r^{\prime}}\right)\,\mathrm{d}\boldsymbol{r^{\prime}}$
(25)

with $k$ and $Z_{0}$ being wavenumber and the free-space impedance,
respectively. The symbol $\overline{\mathbf{G}}$ in (25) denotes dyadic
Green’s function [55] with Cartesian components

${G}_{uv}\left(\boldsymbol{r},\boldsymbol{r^{\prime}}\right)=\left[\delta_{uv}+\frac{1}{k^{2}}\dfrac{\partial^{2}}{\partial
r_{u}\partial
r_{v}}\right]\frac{\mathrm{e}^{-\mathrm{j}k|\boldsymbol{r}-\boldsymbol{r^{\prime}}|}}{4\pi|\boldsymbol{r}-\boldsymbol{r^{\prime}}|},$
(26)

where $\delta_{uv}$ is Kronecker delta and $r_{u}$ denotes the $u$-th
component of radius vector $\boldsymbol{r}$.

The MoM process employs an expansion of unknown surface current density
$\boldsymbol{J}(\boldsymbol{r})$ with a set of basis functions
$\\{\boldsymbol{\psi}_{n}(\boldsymbol{r})\\}$ into

$\boldsymbol{J}(\boldsymbol{r})\approx\displaystyle\sum_{n}I_{n}\boldsymbol{\psi}_{n}(\boldsymbol{r})$
(27)

and Galerkin’s testing procedure [27] to recast (24) into its matrix form

$(\mathbf{Z}_{0}+\mathbf{R}_{\rho})\mathbf{I}=\mathbf{ZI}=\mathbf{V},$ (28)

where $\mathbf{Z}_{0}$ and $\mathbf{R}_{\rho}$ are the vacuum and material
parts of impedance matrix $\mathbf{Z}$, respectively, and $\mathbf{V}$ is the
excitation vector. Explicitly, the operators in (28) are defined element-wise
as

$\displaystyle Z_{0,mn}$
$\displaystyle=\mathrm{j}kZ_{0}\int\limits_{\varOmega}\int\limits_{\varOmega^{\prime}}\boldsymbol{\psi}_{m}(\boldsymbol{r})\cdot\overline{\mathbf{G}}\left(\boldsymbol{r},\boldsymbol{r^{\prime}}\right)\cdot\boldsymbol{\psi}_{n}\left(\boldsymbol{r^{\prime}}\right)\,\mathrm{d}S^{\prime}\,\mathrm{d}S,$
(29) $\displaystyle R_{\rho,mn}$
$\displaystyle=\displaystyle\int\limits_{\varOmega}R_{s}(\boldsymbol{r})\boldsymbol{\psi}_{m}(\boldsymbol{r})\cdot\boldsymbol{\psi}_{n}(\boldsymbol{r})\,\mathrm{d}S,$
(30) $\displaystyle V_{m}$
$\displaystyle=\int\limits_{\varOmega}\boldsymbol{\psi}_{m}(\boldsymbol{r})\cdot\boldsymbol{E}^{\text{i}}\left(\boldsymbol{r}\right)\,\mathrm{d}S.$
(31)

Without loss of generality, we utilize the triangular discretization and Rao-
Wilton-Glisson basis functions [35] which are attached to each inner edge and
its two adjacent triangles. Within our design parametrization, the triangles
have variable surface resistivity $R_{\mathrm{s}}(\rho_{t})$ according to the
interpolation scheme (6). Hence, we can rewrite (30) as

$\mathbf{R}_{\rho}=\sum_{t=1}^{T}R_{\mathrm{s}}(\rho_{t})\mathbf{\Psi}_{\rho}^{t},$
(32)

where $\mathbf{\Psi}_{\rho}^{t}$ is a material element matrix that contains
connectivity properties of a specific triangle $t$. This matrix is written as

$\Psi_{\rho,mn}^{t}=\int\limits_{A_{t}}\boldsymbol{\psi}_{m}\cdot\boldsymbol{\psi}_{n}\,\mathrm{d}A_{t},$
(33)

where the indices $m$ and $n$ uniquely label the basis functions, forming a
block matrix as shown in Fig. 13. The integral in (33) is evaluated
analytically [56].

Figure 13: The example of a particular material element matrix for triangle
$t$ is presented. This triangle is shared by three basis functions
$\\{\boldsymbol{\psi}_{1},\boldsymbol{\psi}_{2},\boldsymbol{\psi}_{3}\\}$,
each of which also occupies one adjacent triangle. These basis functions have
unique labels (given by the numbering of edges) which determine certain
positions with non-zero entries in element matrix $\mathbf{\Psi}_{\rho}^{t}$.
The evaluation of the integral with RWG basis function can be found in [56].

##  Appendix B Influence of the Interpolation scheme on Q-factor

The behavior of the interpolation scheme (6) on the Q-factor at electrical
size $ka=0.59$ is investigated using a model of a loop antenna. Physically
accurate and admissible scenarios are discussed first. Assuming a loop antenna
made of PEC close to resonance, i.e., $Q_{\mathrm{e}}\approx Q_{\mathrm{m}}$,
the black crosses at $\rho=1$ in Fig. 14 represent this scenario. Furthermore,
by considering only the black portion of the structure in Fig. 14 made of PEC,
with the remainder being a vacuum, an electrically short dipole is formed, and
the black crosses indicate the corresponding Q-factor values at $\rho=0$ in
Fig. 14.

The effect of the interpolation scheme (6) can be presented by a continuous
change of surface resistivity $R_{\mathrm{s}}(\rho)$ in the gray part of Fig.
14. Q-factor $Q_{\mathrm{e/m}}$ is evaluated for various values of the lower
bound $R_{\mathrm{s}}(\rho=0)=R_{\mathrm{vac}}$ and upper bound
$R_{\mathrm{s}}(\rho=1)=R_{\mathrm{met}}$ on the attainable resistivities, as
shown in Fig. 14. It can be observed that the physically admissible choice of
boundary resistivities is represented by the green and blue lines in Fig. 14,
which reliably reflect the physical behavior for $\rho=1$ and $\rho=0$. On the
other hand, the combination represented by the red lines does not accurately
reflect the physics at all. The blue and green lines differ mainly in the
intermediate region of $\rho$, where a small change in $\rho$ should produce a
measurable effect on the Q-factor evaluation. The blue dashed line, i.e., a
combination of $[R_{\mathrm{vac}},R_{\mathrm{met}}]=[10^{5},1]\;\Omega$,
delivers the most considerable freedom in the design variable for the
optimizer while minimizing the Q-factor.

Selecting $[R_{\mathrm{vac}},R_{\mathrm{met}}]=[10^{5},1]\;\Omega$ as boundary
resistivity values results in an inaccurate modeling of the physical behavior.
However, based on the authors’ observations while analyzing an optimized
design with proper resistivity values, such as
$[R_{\mathrm{vac}},R_{\mathrm{met}}]\approx[\infty,0]\;\Omega$, there is only
a negligible drop in Q-factor performance. Therefore, it is advantageous to
ease the optimization process at the expense of the small accuracy loss in
Q-factor evaluation. It is important to stress that this claim is not
universal and depends on the interpolation scheme or optimized metric. Hence,
analogous investigations are always recommended.

Figure 14: The influence of interpolation scheme (6) on the electric and
magnetic quality factor $Q_{\mathrm{e/m}}$ is examined at electrical size
$ka=0.59$. The loop antenna, self-resonant at the aforementioned $ka$, is
parametrized by variable $\rho$ representing surface resistivity
$R_{\mathrm{s}}(\rho)$ through (6). The structure is excited by a discrete
delta gap placed in the middle of the left boundary. The antenna is divided
into a fixed part with $\rho=1$ and a designable part with $\rho\in[0,1]$. The
solid lines in the plot show the values of Q-factor $Q_{\mathrm{e/m}}$ for
various combinations of the lower bounds $R_{\mathrm{s}}(0)=R_{\mathrm{vac}}$
as variable $\rho$ is swept, with a fixed value of the upper bound
$R_{\mathrm{s}}(1)=R_{\mathrm{met}}=10^{-2}\;\Omega$. Dashed lines depict the
same sweep but with the upper bound value set to $R_{\mathrm{met}}=1\;\Omega$.

##  Appendix C Q-Factor and its Adjoint Sensitivi ty Analysis

The radiation Q-factor of an antenna can be defined according to [57] as

$Q=\frac{2\omega\max\\{W_{\mathrm{e}},W_{\mathrm{m}}\\}}{P_{\mathrm{rad}}}\equiv\max\\{Q_{\mathrm{e}},Q_{\mathrm{m}}\\},$
(34)

where $\omega$ is angular frequency, $W_{\mathrm{e}}$ and $W_{\mathrm{m}}$ are
cycle mean stored electric and magnetic energy, respectively, and
$P_{\mathrm{rad}}$ is the cycle mean radiated power

$P_{\mathrm{rad}}=\frac{1}{2}\mathbf{I}^{\mathrm{H}}\mathbf{R}_{0}\mathbf{I},$
(35)

where $\mathbf{R}_{0}$ is a real part of the vacuum impedance matrix and
superscript H is the complex conjugate transpose. Stored energies in (34) are
approximated [58], [59] as

$\displaystyle W_{\mathrm{e}}$
$\displaystyle\approx\frac{1}{8}\mathbf{I}^{\mathrm{H}}\left(\frac{\partial\mathbf{X}_{0}}{\partial\omega}-\frac{\mathbf{X}_{0}}{\omega}\right)\mathbf{I}=\frac{1}{4\omega}\mathbf{I}^{\mathrm{H}}\mathbf{X}_{\mathrm{e}}\mathbf{I},$
(36) $\displaystyle W_{\mathrm{m}}$
$\displaystyle\approx\frac{1}{8}\mathbf{I}^{\mathrm{H}}\left(\frac{\partial\mathbf{X}_{0}}{\partial\omega}+\frac{\mathbf{X}_{0}}{\omega}\right)\mathbf{I}=\frac{1}{4\omega}\mathbf{I}^{\mathrm{H}}\mathbf{X}_{\mathrm{m}}\mathbf{I},$
(37)

$\mathbf{X}_{0}$ is an imaginary part of the vacuum impedance matrix, and
$\mathbf{X}_{0}=\mathbf{X}_{\mathrm{m}}-\mathbf{X}_{\mathrm{e}}$. These matrix
operators are obtained from method of moments codes [47]. The definition (34)
is interpreted as the upper bound on the Q-factor of an antenna system that is
tuned to the resonance by the addition of a lossless reactive element.

The sensitivities of $Q_{\mathrm{e}}$ and $Q_{\mathrm{m}}$ with respect to
varible $\rho_{t}$ can be computed by adjoint sensitivity analysis [14] as

$\frac{\,\mathrm{d}Q_{\mathrm{e/m}}}{\,\mathrm{d}\rho_{t}}=\frac{\partial
Q_{\mathrm{e/m}}}{\partial\rho_{t}}+2\mathrm{Re}\left\\{\boldsymbol{\lambda}_{\mathrm{e/m}}^{\mathrm{T}}\frac{\partial\mathbf{R}_{\rho}}{\partial\rho_{t}}\mathbf{I}\right\\},$
(38)

where $\boldsymbol{\lambda}_{\mathrm{e/m}}$ is the solution to corresponding
adjoint equations

$\mathbf{Z}^{\mathrm{T}}\boldsymbol{\lambda}_{\mathrm{e/m}}=-\left(\frac{\partial
Q_{\mathrm{e/m}}}{\partial\mathbf{I}}\right)^{\mathrm{T}},$ (39)

where the right-hand side reads

$\frac{\partial
Q_{\mathrm{e/m}}}{\partial\mathbf{I}}=\frac{1}{2P_{\mathrm{rad}}}\mathbf{I}^{\mathrm{H}}\left(\mathbf{X}_{\mathrm{e/m}}-Q_{\mathrm{e/m}}\mathbf{R}_{0}\right).$
(40)

The Q-factor is not a function of $\rho_{t}$, and the first term in (38)
vanishes.

##  Appendix D Resilience towards initial density distribution

The Q-factor minimization by the density-based topology optimization on a
rectangular design region from Section V-A is considered here with the same
settings but different initial density distribution, namely, uniform
distribution with $\rho_{t}=0.1$, random distribution, and distribution with
already well-separated charge distribution (low Q-factor value), see the
bottom of Fig. 15. The convergence curves are reported at the top of Fig. 15.
The method delivered optimized designs, see bottom of Fig. 15, with similar
values of Q-factor $Q/Q_{\mathrm{lb}}^{\mathrm{TM}}=1.06$, i.e., the nature of
Q-factor and usage of spatial filtering make the proposed method less
dependent on the initial seed. The main influence of the initial seed lies in
the convergence sharpness at the early iterations. This study might
nonetheless be not applicable to EM metrics different from the Q-factor.

Figure 15: (Top) Convergence history of Q-factors
$Q=\mathrm{max}\\{Q_{\mathrm{e}},Q_{\mathrm{m}}\\}$ normalized to the
fundamental bound $Q_{\mathrm{lb}}^{\mathrm{TM}}$ [45] for three different
initial density distribution. Occasional jumps result from the projection
filter update. (Bottom) Snapshots of an initial and optimized density field
distribution are included below in the convergence graph. The outline of each
design is color-matched with the curve in the convergence plot.

##  Appendix E Parametrization of Spherical Curves

A spherical helix [60] is a curve maintaining a constant angle with the
equatorial plane and is described by the number of turns $M$ which connect the
poles. A parametrization of a spherical helix is given by

$\boldsymbol{r}(t)=R\Big{[}f(t)\cos(2\pi Mt),f(t)\sin(2\pi Mt),2t-1\Big{]},$
(41)

where $f(t)=\sqrt{1-4(t-0.5)^{2}}$, $R$ is the radius of the sphere and where
different components of (41) refer to Cartesian coordinates.

Loxodrome (Rhumb line) [54] is a curve that crosses all meridians at the same
angle. Its parametrization in Cartesian coordinates reads

$\boldsymbol{r}(t)=R\left[\dfrac{\cos(t)}{\cosh(\gamma
t)},\dfrac{\sin(t)}{\cosh(\gamma t)},\tanh(\gamma t)\right],$ (42)

where $\gamma$ is the slope of the curve. Alternative parametrization in
spherical coordinates is given in [31]. Denoting
$\boldsymbol{\tau}=\mathrm{d}\boldsymbol{r}/\mathrm{d}t$ a tangent vector
field along the curve, it can be shown that

$\left|\boldsymbol{\tau}\right|=\dfrac{\sqrt{1+\gamma^{2}}}{\cosh\left(\gamma
t\right)},$ (43)

and that spherical vector components of the unit vector along the curve read

$\dfrac{\tau_{r}}{\left|\boldsymbol{\tau}\right|}=0,\,\dfrac{\tau_{\theta}}{\left|\boldsymbol{\tau}\right|}=-\dfrac{\gamma}{\sqrt{1+\gamma^{2}}},\,\dfrac{\tau_{\varphi}}{\left|\boldsymbol{\tau}\right|}=\dfrac{1}{\sqrt{1+\gamma^{2}}}.$
(44)

This is the vector orientation of an addition of $\mathrm{TM}_{10}$ and
$\mathrm{TE}_{10}$ spherical vector waves, i.e.., the correct orientation to
minimize Q-factor on a spherical shell [31], [53].

## References

  * [1] C. A. Balanis, _Antenna Theory Analysis and Design_ , 3rd ed. Wiley, 2005.
  * [2] K. Fujimoto and H. Morishita, _Modern Small Antennas_. Cambridge, Great Britain: Cambridge University Press, 2013\.
  * [3] J. L. Volakis, C. Chen, and K. Fujimoto, _Small Antennas: Miniaturization Techniques & Applications_. McGraw-Hill, 2010.
  * [4] S. Koziel and S. Ogurtsov, _Antenna design by simulation-driven optimization_. Springer, 2014.
  * [5] S. Koziel and A. Pietrenko-Dabrowska, “Robust parameter tuning of antenna structures by means of design specification adaptation,” _IEEE Transactions on Antennas and Propagation_ , vol. 69, no. 12, pp. 8790–8798, 2021\.
  * [6] J. Robinson and Y. Rahmat-Samii, “Particle swarm optimization in electromagnetics,” _IEEE Trans. Antennas Propag._ , vol. 52, no. 2, pp. 397–407, Feb. 2004.
  * [7] S. Koziel and A. Pietrenko-Dabrowska, “Improved-efficacy em-driven optimization of antenna structures using adaptive design specifications and variable-resolution models,” _IEEE Transactions on Antennas and Propagation_ , 2023.
  * [8] M. P. Bendsoe and O. Sigmund, _Topology Optimization_ , 2nd ed. Berlin, Germany: Springer, 2004.
  * [9] Y. Rahmat-Samii and E. Michielssen, Eds., _Electromagnetic Optimization by Genetic Algorithms_. Wiley, 1999.
  * [10] M. Capek, M. Gustafsson, L. Jelinek, and P. Kadlec, “Optimal Inverse Design Based on Memetic Algorithms – Part 1: Theory and Implementation,” eprint arXiv: 2110.08044v2. [Online]. Available: https://arxiv.org/abs/2110.08044v2
  * [11] ——, “Optimal Inverse Design Based on Memetic Algorithms – Part 2: Examples and Properties,” eprint arXiv: 2110.13460v2. [Online]. Available: https://arxiv.org/abs/2110.13460v2
  * [12] G. L. Nemhauser and L. A. Wolsey, _Integer an Combinatorial Optimization_. John Wiley & Sons, 1999\.
  * [13] O. Sigmund, “On the usefulness of non-gradient approaches in topology optimization,” _Struct. Multidisc. Optim._ , vol. 43, pp. 589–596, 2011\.
  * [14] D. A. Tortorelli and P. Michaleris, “Design sensitivity analysis: overview and review,” _Inverse problems in Engineering_ , vol. 1, no. 1, pp. 71–105, 1994\.
  * [15] 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.
  * [16] M. P. Bendsøe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” _Computer methods in applied mechanics and engineering_ , vol. 71, no. 2, pp. 197–224, 1988.
  * [17] O. Sigmund and K. Hougaard, “Geometric properties of optimal photonic crystals,” _Physical review letters_ , vol. 100, no. 15, p. 153904, 2008\.
  * [18] J. S. Jensen and O. Sigmund, “Systematic design of photonic crystal structures using topology optimization: Low-loss waveguide bends,” _Applied physics letters_ , vol. 84, no. 12, pp. 2022–2024, 2004.
  * [19] T. Nomura, K. Sato, K. Taguchi, T. Kashiwa, and S. Nishiwaki, “Structural topology optimization for the design of broadband dielectric resonator antennas using the finite difference time domain technique,” _International Journal for Numerical Methods in Engineering_ , vol. 71, no. 11, pp. 1261–1296, 2007.
  * [20] A. Erentok and O. Sigmund, “Topology optimization of sub-wavelength antennas,” _IEEE Trans. Antennas Propag._ , vol. 59, no. 1, pp. 58–69, 2011\.
  * [21] N. Aage and V. Egede Johansen, “Topology optimization of microwave waveguide filters,” _International Journal for Numerical Methods in Engineering_ , vol. 112, no. 3, pp. 283–300, 2017.
  * [22] E. Hassan, E. Wadbro, and M. Berggren, “Topology optimization of metallic antennas,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 5, pp. 2488–2500, May 2014.
  * [23] Q. Wang, R. Gao, and S. Liu, “A novel parameterization method for the topology optimization of metallic antenna design,” _Acta Mechanica Sinica_ , vol. 33, no. 6, pp. 1040–1050, 2017.
  * [24] 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, 2001.
  * [25] M. Cismasu and M. Gustafsson, “Antenna bandwidth optimization by genetic algorithms with single frequency simulation,” in _Proceedings of the 7th European Conference on Antennas and Propagation (EUCAP)_ , Apr. 2013, pp. 2781–2782.
  * [26] 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\.
  * [27] R. F. Harrington, _Field Computation by Moment Methods_. Piscataway, New Jersey, United States: Wiley – IEEE Press, 1993.
  * [28] 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, 2011.
  * [29] 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.
  * [30] 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,” vol. 61, no. 4, pp. 19–29, Aug. 2019.
  * [31] O. S. Kim, “Minimum Q electrically small antennas,” _IEEE Trans. Antennas Propag._ , vol. 60, no. 8, pp. 3551–3558, Aug. 2012.
  * [32] 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.
  * [33] “Topology optimization for minimal Q-factor.” [Online]. Available: {https://github.com/tucekjon/Topology_optimization_minQ}
  * [34] W. C. Gibson, _The Method of Moments in Electromagnetics_ , 2nd ed. Chapman and Hall/CRC, 2014.
  * [35] 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.
  * [36] M. P. Bendsøe and O. Sigmund, “Material interpolation schemes in topology optimization,” _Archive of applied mechanics_ , vol. 69, no. 9, pp. 635–654, 1999.
  * [37] G. Kiziltas, D. Psychoudakis, J. L. Volakis, and N. Kikuchi, “Topology design optimization of dielectric substrates for bandwidth improvement of a patch antenna,” _IEEE Transactions on Antennas and Propagation_ , vol. 51, no. 10, pp. 2732–2743, 2003.
  * [38] R. L. Fante, “Quality factor of general ideal antennas,” _IEEE Trans. Antennas Propag._ , vol. 17, no. 2, pp. 151–157, 1969.
  * [39] L. J. Chu, “Physical limitations of omni-directional antennas,” _J. Appl. Phys._ , vol. 19, pp. 1163–1175, Dec. 1948.
  * [40] A. Diaz and O. Sigmund, “Checkerboard patterns in layout optimization,” _Structural optimization_ , vol. 10, no. 1, pp. 40–45, 1995.
  * [41] R. B. Haber, C. S. Jog, and M. P. Bendsøe, “A new approach to variable-topology shape design using a constraint on perimeter,” _Structural optimization_ , vol. 11, no. 1, pp. 1–12, 1996.
  * [42] O. Sigmund, “On the design of compliant mechanisms using topology optimization,” _Journal of Structural Mechanics_ , vol. 25, no. 4, pp. 493–524, 1997.
  * [43] B. Bourdin, “Filters in topology optimization,” _International journal for numerical methods in engineering_ , vol. 50, no. 9, pp. 2143–2158, 2001.
  * [44] 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, 1987.
  * [45] 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. [Online]. Available: http://www.e-fermat.org/articles/gustafsson-art-2016-vol15-may-jun-005/
  * [46] (2022) The Matlab. The MathWorks. [Online]. Available: www.mathworks.com
  * [47] (2022) Antenna Toolbox for MATLAB (AToM). Czech Technical University in Prague. www.antennatoolbox.com. [Online]. Available: www.antennatoolbox.com
  * [48] M. Gustafsson, C. Sohl, and G. Kristensson, “Illustrations of new physical bounds on linearly polarized antennas,” _IEEE Trans. Antennas Propag._ , vol. 57, no. 5, pp. 1319–1327, May 2009.
  * [49] A. Clausen and E. Andreassen, “On filter boundary conditions in topology optimization,” _Structural and Multidisciplinary Optimization_ , vol. 56, pp. 1147–1155, 2017.
  * [50] M. Gustafsson, M. Cismasu, and B. L. G. Jonsson, “Physical bounds and optimal currents on antennas,” LUTEDX/(TEAT-7210)/1-22/, Tech. Rep., 2011.
  * [51] A. R. Diaz and O. Sigmund, “A topology optimization method for design of negative permeability metamaterials,” _Structural and Multidisciplinary Optimization_ , vol. 41, no. 2, pp. 163–177, 2010.
  * [52] N. Aage, N. Mortensen, and O. Sigmund, “Topology optimization of metallic devices for microwave applications,” _International Journal for numerical methods in engineering_ , vol. 83, no. 2, pp. 228–248, 2010.
  * [53] 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.
  * [54] J. Alexander, “Loxodromes: A rhumb way to go,” _Mathematics magazine_ , vol. 77, no. 5, pp. 349–356, 2004.
  * [55] W. C. Chew, _Waves and Fields in Inhomogeneous Media_. IEEE Press, 1995.
  * [56] L. Jelinek and M. Capek, “Optimal currents on arbitrarily shaped surfaces,” _IEEE Trans. Antennas Propag._ , vol. 65, no. 1, pp. 329–341, Jan. 2017.
  * [57] R. E. Collin and S. Rothschild, “Evaluation of antenna Q,” _IEEE Trans. Antennas Propag._ , vol. 12, no. 1, pp. 23–27, Jan. 1964.
  * [58] G. A. E. Vandenbosch, “Reactive energies, impedance, and Q factor of radiating structures,” _IEEE Trans. Antennas Propag._ , vol. 58, no. 4, pp. 1112–1127, Apr. 2010.
  * [59] M. Gustafsson and S. Nordebo, “Optimal antenna currents for Q, superdirectivity, and radiation patterns using convex optimization,” _IEEE Trans. Antennas Propag._ , vol. 61, no. 3, pp. 1109–1118, Mar. 2013\.
  * [60] P. D. Scofield, “Curves of constant precession,” _The American mathematical monthly_ , vol. 102, no. 6, pp. 531–537, 1995.

|  Jonas Tucek received the M.Sc. degree in electrical engineering from Czech
Technical University, Prague, Czech Republic, in 2021, where he is currently
pursuing the Ph.D. degree in the area of topology optimization.
---|---
|  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, antenna design, numerical techniques, and
optimization. He authored or co-authored over 130 journal and conference
papers. Dr. Capek is the Associate Editor of IET Microwaves, Antennas &
Propagation. He was a regional delegate of EurAAP between 2015 and 2020 and an
associate editor of Radioengineering between 2015 and 2018.
---|---
|  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.
---|---
|  Ole Sigmund received the Ph.D. degree and Habilitation in 1994 and 2001,
respectively, and has had research positions with the University of Essen and
Princeton University. He is a professor in the Department of Mechanical
Engineering, Technical University of Denmark (DTU). He is a member of the
Danish Academy of Technical Sciences and the Royal Academy of Science and
Letters, Denmark and is the former elected President (2011-15, now EC member)
of ISSMO (International Society of Structural and Multidisciplinary
Optimization). His research interests include theoretical extensions and
applications of topology optimization methods to mechanics, and multiphysics
problems.
---|---