# An Improved Algorithm for Coarse-Graining Cellular Automata

Yerim Song  Canyon Crest Academy, San Diego, CA, USA<EMAIL_ADDRESS>Joshua A. Grochow  Department of Computer Science, University of Colorado
Boulder, Boulder, CO, USA  Department of Mathematics, University of Colorado
Boulder, Boulder, CO, USA<EMAIL_ADDRESS>
###### Abstract

In studying the predictability of emergent phenomena in complex systems,
Israeli & Goldenfeld (_Phys. Rev. Lett._ , 2004; _Phys. Rev. E_ , 2006) showed
how to coarse-grain (elementary) cellular automata (CA). Their algorithm for
finding coarse-grainings of supercell size $N$ took doubly-exponential
$2^{2^{N}}$-time, and thus only allowed them to explore supercell sizes $N\leq
4$. Here we introduce a new, more efficient algorithm for finding coarse-
grainings between any two given CA that allows us to systematically explore
all elementary CA with supercell sizes up to $N=7$, and to explore individual
examples of even larger supercell size. Our algorithm is based on a
backtracking search, similar to the DPLL algorithm with unit propagation for
the $\mathsf{NP}$-complete problem of Boolean Satisfiability.

##  I Introduction

Cellular automata (CA) are a model of dynamical systems that are discrete in
both space and time. Since their development by Ulam and von Neumann [1] in
the 1940s, CA have been applied to many subjects including biology [2],
physics [3], and computer science [4, 5].

Coarse-graining is one method of reducing the complexity of a system, while
hopefully retaining enough structure to reveal insights about it. Rather than
focusing on small, specific details, we “zoom out” and look at a larger
picture. Israeli & Goldenfeld [6, 7] introduced a method of coarse-graining
cellular automata, and applied it systematically to the _elementary_ cellular
automata.

Elementary cellular automata were introduced and studied extensively by
Wolfram [6]: they are automata with 2 states, on a 1-dimensional grid of
cells, with only nearest-neighbor interactions. There are exactly 256 distinct
elementary CAs, and Wolfram proposed dividing them into four classes based on
their long-term behavior: stabilizing to a homogeneous state, ultimately
becoming periodic, or maintaining random or complex-looking behavior
indefinitely. The latter classes are hypothesized to be “computationally
irreducible” and difficult or impossible to predict [8, 9, 10]. However
through the process of coarse-graining, even CA that are originally in these
complex classes can coarse-grain to CA in the simpler classes [6, 7].

To calculate these coarse-grainings, Israeli and Goldenfeld used a brute force
algorithm, which required $2^{2^{N}}$ steps to search for coarse-grainings of
supercell size $N$. The lengthy time required for the computations limited
their study to coarse-grainings of supercell size only 4. They reported that
there were 16 elementary CA that they could not coarse-grain at all, and for
the other CA it was unknown whether they had found all the available coarse-
grainings.

In this paper we develop a new algorithm, using a backtracking search similar
to the DPLL algorithm for Boolean Satisfiability [11, 12]. This backtracking
search allows us to prune branches from the search tree early, based on
coarse-graining constraints. Furthermore, dynamically ordering the variables
allows us to speed up the process considerably in practice.

Using our new algorithm, we were able to systematically find all coarse-
grainings between elementary CA up to supercell size 7 on a commodity laptop.
Note that $2^{2^{7}}\approx 3.4\times 10^{38}$, so even if each step could be
done in 1ms the brute force algorithm of [6, 7] would take $\approx 10^{27}$
years to handle supercell size 7. We were also able to search for specific
coarse-grainings of larger supercell size. We find 56 new coarse-grainings at
these larger supercell sizes (26 up to symmetry), though we leave open the
mathematical question of how to prove when all nontrivial coarse-grainings
have been found. Interestingly, all of these new coarse-grainings were found
at $N=5$ or $6$; there were no new pairs of CA $A,B$ such that $A$ coarse-
grained to $B$ at supercell size $7$.

###  I.1 Related Work

Magiera & Dzwinel [13, 14] also present an improved algorithm for coarse-
graining CA, and get up to supercell size 7. However, despite the similar
name, they are solving a related but different problem than the one we solve.
Namely, they solve the problem:

> IsCoarseGrainable
>
> Input: A CA $A$ and supercell size $N$
>
> Output: A nontrivial coarse-graining of $A$ with supercell size $N$, or
> $\bot$ if none exists.

In contrast, we solve the following problem, which enables us to fill out Fig.
5 (a CA analogue of a renormalization group flow diagram):

> CoarseGraining
>
> Input: Two CAs $A,B$ and supercell size $N$
>
> Output: A nontrivial coarse-graining from $A$ to $B$ with supercell size
> $N$, or $\bot$ if none exists.

Their algorithm begins by constructing $A^{N}$, with $k^{N}$ states per cell
if $A$ had $k$ states per cell, and then tries to collapse those new states
together as much as possible, resulting in another CA (but one which is not
specified ahead of time in the input). In contrast, our algorithm tries to
find the coarse-graining map from $A^{N}$ to the given CA $B$. The fact that
$B$ is given in the input allows us to use other algorithmic tactics that are
less easy to take advantage of in their setting. In contrast, their algorithm
is aimed at deciding which CA can be coarse-grained at all.

We note that one could directly reduce to a Boolean Satisfiability problem and
use an off-the-shelf SAT solver. This did not seem very promising to us, as
supercell size $N$ would result in a SAT instance with $2^{N}$ variables and
$2^{3N}$ constraints, where each constraint would consist of several CNF
clauses. See Remark 1 below for more details.

However, given that our approach is similar to approaches to the
$\mathsf{NP}$-complete problem of Boolean Satisfiability, it is natural to
wonder about comparing the two. Because the coarse-graining problem has
exponentially many variables—one variable for the projection of each possible
$N$-tuple of states—it is easy to see that CoarseGraining is in the complexity
class $\mathsf{NEXP}$. It would be interesting to know whether it is
$\mathsf{NEXP}$-complete, as that would suggest that some amount of brute
force search is inevitable, assuming the widely believed complexity conjecture
that $\mathsf{EXP}\neq\mathsf{NEXP}$.

We note that while we present our algorithm and results in the context of
_elementary_ CA (2-state, 1D, nearest neighbor), it is trivial to adapt it to
arbitrary CA.

##  II Background on Cellular Automata and Coarse-Graining

A cellular automaton (CA) $A$ has cells at the nodes of some (usually
infinite, regular) graph, each cell has a state from a finite set $S_{A}$, and
an update rule $f_{A}(a_{0};a_{1},\dotsc,a_{k})$, which takes in the state
$a_{0}$ of a cell and the states $a_{1},\dotsc,a_{k}$ of its neighbors, and
outputs the state of the cell for the next time step. In this paper we
consider only _elementary_ CA, which are CA on a 1D lattice, with nearest-
neighbor interactions, and 2 states, $S_{A}=\\{0,1\\}$. Thus each rule $f_{A}$
takes in only three states: that of a cell, its left neighbor, and its right
neighbor. For geometric convenience, rather than putting the cell’s state
first, we put them in geometric order, so that, e.g.,
$f_{A}(a_{1},a_{2},a_{3})$ is the next state of the cell in position 2.

As introduced by Israeli & Goldenfeld [6, 7], a coarse-graining from one CA
$A$ to another CA $B$ with supercell size $N$ is a map $P\colon S_{A}^{N}\to
S_{B}$ such that, given any string states, if we first run $A$ for $N$ steps
and then apply $P$, the result is the same as first applying $P$, and then
running $B$ for a single step. That is, the coarse-graining scales both the
cell size and the speed of the CA by $N$.

For 1D CA with nearest-neighbor interactions, this is captured precisely by
the fundamental coarse-graining equation

$P(f_{A}^{N}(x_{1},x_{2},x_{3}))=f_{B}(P(x_{1}),P(x_{2}),P(x_{3}))\quad\forall
x_{i}\in S_{A}^{N}$ (1)

Given an input $(x_{1},x_{2},x_{3})\in S_{A}^{3N}$, we will refer to the
result on the left-hand side here throughout the paper as
$res_{1}=res_{1}(x_{1},x_{2},x_{3})$, and similarly the result on the right-
hand side as $res_{2}$. We will sometimes abbreviate Eq. (1) to:

$PA^{N}=BP.$

Given CA $A,B$ and supercell size $N$, our goal is thus to find a projection
$P:S_{A}^{N}\to S_{B}$ such that

$res_{1}=res_{2}\quad\forall(x_{1},x_{2},x_{3})\in(S_{A}^{N})^{3}.$

These will be the fundamental constraints that our backtracking algorithm
explores and exploits.

Following [6, 7], we consider a coarse-graining “trivial” if $P$ is a constant
function, that is, if it maps all inputs to $0$ (resp., all inputs to $1$). We
ignore these trivial coarse-grainings by fiat.

It is difficult to discover which projection functions $P$ are valid coarse-
grainings at window size $N$, because there are $2^{2^{N}}-2$ possibilities
for $P$ (excluding the two trivial possibilities), and for each such $P$ there
are $2^{3N}$ constraints to check. There seems to be no apparent pattern as to
which projections will satisfy all of constraints, and thus be a coarse-
graining. (The lack of such patterns could be attested to theoretically by
answering affirmatively our question above about
$\mathsf{NEXP}$-completeness.) In the absence of such patterns, one might be
resigned to trying every possible projection, and indeed this is essentially
the approach taken in [6, 7]. In this paper we show that a backtracking search
with a few easy-to-implement heuristics can do significantly better in
practice, while still ensuring our algorithm is _complete_ in the sense that
every possible coarse-graining will be found.

###  II.1 Symmetries

In our algorithm, we will determine whether there are coarse-grainings from
$A$ to $B$ for all pairs of elementary CA $A,B$, for supercell sizes up to
$N=7$. Once we have found a coarse-graining from $A$ to $B$, we do not
consider it further in our experiments; e. g., if we find a coarse-graining
from $A$ to $B$ at supercell size $3$, we do not search for coarse-grainings
from $A$ to $B$ with larger supercell sizes, though our algorithm could be
used to do so for further exploratory purposes.

Additionally, as noted in [6, 7], we may eliminate CA that are “isomorphic” to
one another under the symmetries that swap left and right, that swap 0 and 1,
or that swap both left-right and 0-1. If we use $\sigma_{LR}$ to denote the
left-right swap, $\sigma_{01}$ to denote the $0$-$1$ swap, and $A\to B$ to
denote a coarse-graining, then we have:

$\displaystyle A\to B\Longleftrightarrow$ $\displaystyle\sigma(A)\to\sigma(B)$
$\displaystyle\forall\sigma\in\\{\sigma_{LR},\sigma_{01},\sigma_{LR}\circ\sigma_{01}=\sigma_{01}\circ\sigma_{LR}\\}.$

For example, suppose $A$ is rule 24 and $B$ is rule 240. We have
$\sigma_{01}(24)=231,\sigma_{LR}(24)=66,\sigma_{01,LR}(24)=189$, and
$\sigma_{01}(240)=240,\sigma_{LR}(240)=\sigma_{LR,01}(240)=170$. So among the
possibilities

$24\to 240\qquad 231\to 240\qquad 66\to 170\qquad 189\to 170,$

we need only check one, rather than all four. By eliminating these symmetric
cases, the number of comparisons we must make between pairs of elementary CA
is cut down to about 1/3 of the original total of $256^{2}$; it is not 1/4
because of the presence of some rules that are mapped to themselves by some of
the symmetries (such as $\sigma_{01}(240)=240$ in the preceding example).

##  III Backtracking Algorithm

The new approach we developed to more efficiently compute coarse-grainings
involves considering a backtracking tree search, allowing us to to selectively
cut tree branches earlier in the computation without losing information. The
basic idea is that, given cellular automata $A,B$ and supercell size $N$, we
seek to find values for $P(x)$ for all $x\in S_{A}^{N}$, that satisfy (1). For
each such $x$, we branch on trying either $P(x)=0$ or $P(x)=1$ (when working
with elementary CA). As we assign values to $P(x)$ for some of the $x$’s, we
can begin to see whether (1) is satisfied or violated for different inputs
$(x_{1},x_{2},x_{3})\in S_{A}^{3N}$. In some cases, this lets us undo an
assignment for $P(x)$, backtrack, and try the other value.

Figure 1 shows an example of this search tree for supercell size $N=2$.
Generally, for supercell size $N$, the tree will have $2^{N}$ levels and
$2^{2^{N}}$ leaves, though many of these will be pruned before they are
reached. Note that although it is included in the tree diagram, we exclude the
all-0s and all-1s projections in our calculations, as these are considered
“trivial” coarse-grainings, as in Israeli & Goldenfeld. These projections are
represented in the first and last leaves in the diagram. If we did not exclude
them, then almost all rules could be coarse grained by these trivial
projections.

Figure 1:  The search tree above displays all 16 possible projections with
supercell size $N=2$. Based on the projection input on the left, the output is
displayed on the tree. In this diagram, level 1 represents $P(00)$, level 2
represent $P(01)$, and so on. The boxed sequence demonstrates that the fourth
leaf corresponds to the projection $0011$, that is,
$P(00)=0,P(01)=0,P(10)=1,P(11)=1$.

###  III.1 Tree pruning

The first way we reduce the number of computations is to prune branches at
earlier levels of the tree, so that computations in the remaining levels can
be skipped. By allowing several projections to be skipped early on, we manage
to significantly reduce the computation time to find or rule out a coarse-
graining.

To determine whether a branch can be pruned the algorithms checks (1) on
inputs $(x_{1},x_{2},x_{3})\in(S_{A}^{N})^{3}$ that have been determined by
the choices made so far. Given a partial projection $P$, defined only on a
subset of $S_{A}^{N}$, we say that an input $(x_{1},x_{2},x_{3})$ has a
decided projection if $P(x_{1}),P(x_{2}),P(x_{3})$, and
$P(f_{A}^{N}(x_{1},x_{2},x_{3}))$ have all been assigned already. We begin at
level 1, where the first projection will be set to either 0 or 1. Meanwhile,
the rest of the projections will be temporarily undecided. Using all possible
inputs with decided projections, we check if $PA^{N}=BP$ is satisfied before
proceeding to the next level. Recall that we use $res_{1}$ to denote the
output of $PA^{N}$ (when defined) and $res_{2}$ the output of $BP$. If we
discover a decided input for which $res_{1}\neq res_{2}$, then the most recent
branch leading to the current projection can be safely cut off, as there are
no valid projections that extend it (i.e., in its subtree).

######  Example 1.

To better illustrate the method, we walk through an example of finding a
coarse-graining from rule 196 to rule 192 with supercell size $N=2$; this
example is illustrated in Figure 2. To begin, the tree diagram will start at
level 1 and branch on $P(00)=0$. We only check inputs whose projections are
defined, so the only available input to test at this point is $x=000000$. When
rule 196 is applied to this pattern $N=2$ times, the result is $00$, so we may
indeed check the result at $x$. $P(00)=0$, so $res_{1}=0$. To calculate
$res_{2}$, we apply $P$ to the input $000000$ three times, resulting in
$P(00)P(00)P(00)=000$, and then use rule 192 on $000$ to get $res_{2}=0$.
Because $res_{1}=res_{2}$ at this point, we cannot cut this branch, and
continue to the next level in the tree.

Suppose that in the next level we branch on $P(01)=1$ (and we still have
assigned $P(00)=0$). As $x=000100$ has its projection defined, we may attempt
to validate our choices on that input. To get $res_{1}$ we first apply rule
196 twice to get an outcome of $01$ and then $P$ to the result to get
$PA^{N}(000100)=P(01)=1$. To compute $res_{2}$ we apply rule 192 to
$P(00)P(01)P(00)=010$, resulting in $res_{2}=0$. Since $res_{1}\neq res_{2}$,
we can cut off the whole branch with $P(00)=0,P(01)=1$, avoiding further
computations.

######  Remark 1.

This is one of the places our approaches differs from simply reducing to a
Boolean Satisfiability (SAT) problem and using an off-the-shelf SAT solver. In
our approach, we get to choose “on the fly” which inputs $(x_{1},x_{2},x_{3})$
to test for $res_{1}=res_{2}$. In contrast, a reduction to a SAT solver would
require writing down _all_ the constraints for every possible input $x$ right
from the beginning (of which there are $8^{N}$), before any branches have been
made.

Figure 2:  Tree pruning. The search tree for attempting to coarse grain from
rule 196 to rule 192 at supercell size $N=2$ (Example 1). After branching on
$P(00)=0$ and $P(01)=1$, one of the equations (1) is violated, and the search
backs up to try the branch $P(01)=0$ instead. The ‘X’ represents the tree
being pruned at level 2; the box represents all the computations that were
skipped due to this cut.

###  III.2 Forced values

Our second method to reduce computation time is to leverage forced values of
the projection. Namely, once we have branched on assigning $P(x)$ for some
values of $x\in S_{A}^{N}$, we may be forced to assign $P(x^{\prime})$ for
some new, previously unassigned $x^{\prime}$ in order to satisfy (1). (This is
similar to unit propagation in Boolean Satisfiability solvers.) When this
occurs later in the tree, it lets us eliminate duplicated computations; when
it occurs earlier in the tree, it lets us avoid large swaths of the tree.

In this method, instead of simply checking whether $res_{1}=res_{2}$ on inputs
$(x_{1},x_{2},x_{3})$ for which $P(x_{1}),P(x_{2}),P(x_{3})$, _and_
$P(f_{A}^{N}(x_{1},x_{2},x_{3}))$ have all been assigned, we make an effort to
match $res_{1}$ and $res_{2}$ even when these have not yet all been assigned.
As in the method and example above, we begin with an input
$(x_{1},x_{2},x_{3})$ such that $P(x_{1}),P(x_{2})$, and $P(x_{3})$ are all
defined. We can then calculate $res_{2}=f_{B}(P(x_{1}),P(x_{2}),P(x_{3}))$.
Let $x_{4}=f_{A}^{N}(x_{1},x_{2},x_{3})$; note that $res_{1}=P(x_{4})$. If
$P(x_{4})$ is not yet assigned, then the value of $P(x_{4})$ is forced by our
earlier choices to be equal to $res_{2}$.

######  Example 2.

We illustrate the process of forcing on a coarse-graining from rule 2 to rule
4 with $N=2$; see Figure 3. In this example, we will assign variables in
“backwards” order, beginning with $P(11)$ and ending with $P(00)$. We start by
branching on $P(11)=0$. The only input for which $P(x_{1}),P(x_{2}),P(x_{3})$
are all determined is $x=111111$. To find $res_{1}$ we apply rule 2 $N$ times
and get the result of $f_{A}^{N}(x)=00$. However, $P(00)$ won’t be determined
until the 4th level of the tree, and we are only at level 1, so $P(00)$ is
currently unassigned. We compute $res_{2}$ and get a result of 0. To ensure
$res_{1}=res_{2}$, $res_{1}$ must also be 0, so we record $P(00)=0$. Then the
algorithm continues by branching on the value of $P(10)$, then $P(01)$, and
when it gets to the last level, it remembers that the value of $P(00)$ was
already forced to $0$ and does not branch further; at that point it just
checks whether the choices it has made satisfy (1) for all inputs.

Figure 3:  Forced values. This search tree shows the coarse-graining from rule
2 to rule 4 at supercell size $N=2$. While calculating if we can continue
after branching on $P(11)=0$, we find that the projection at level 4, $P(00)$
is forced to be 0. The algorithm thus records that $P(00)=0$ at this point,
but then continues down the search tree as previously.

###  III.3 Dynamic variable ordering

Our third method to reduce computation time is to take further advantage of
the forced values. Namely, when a value is forced, we can be more efficient by
re-ordering the projection inputs dynamically (i. e., at the time). That is,
instead of the levels of the tree being labeled in a static order such as 00,
01, 10, 11 (“forwards”), or 11, 10, 01, 00 (“backwards”), we can instead make
the ordering dynamic as the algorithm proceeds. By considering a forced
variable as soon as it is forced, the algorithm is able to make further
inferences at the time, performing even more pruning and perhaps finding even
more forced variables. This is analogous to unit propagation in SAT solvers.

######  Example 3 (Continuation of Example 2).

See Figure 4. As in Example 2, suppose we branch on $P(11)=0$ and discover
that $P(00)=0$ is forced. Rather than waiting to see $P(00)$ at level 4 as
“originally planned” (in the static, backwards order), the algorithm decides
to make $P(00)$ the immediate next (second) level of the tree. The algorithm
now continues as before. But while doing this, we see that $P(10)$ is also
forced, so level 3 gets dynamically set to branch on $P(10)$. Dynamic ordering
allows cuts to occur earlier in the search tree, significantly reducing the
computation time.

Figure 4:  The dynamic ordering process. See Example 3. (Above) In this
example, the backwards algorithm is used, so it begins by branching on
$P(11)$, while the ordering of the remaining variables is not yet decided
(indicated by the “?”s). (Below) After branching on $P(11)$, the value of
$P(00)$ is forced, so level 2 is the assignments to 00. The process continues
in this way until all the variables are assigned either by branching or being
forced.

###  III.4 Implementation variants

Our dynamic ordering selects the next projection as one that is forced, if
such a projection exists. We still have to explain what inputs $x$ to test
$res_{1}\stackrel{{\scriptstyle?}}{{=}}res_{2}$ on (when there is a choice,
which occurs frequently further along in the tree), and what projection to
choose next when none are forced, including the first projection to branch on.

We implemented four specific strategies for these choices that we test
experimentally:

  * $F_{1}$

In this strategy, the default ordering is “forward,” moving from $P(00\dotsb
0)$ to $P(11\dotsb 1)$, treating the input strings as binary representations
of integers in the usual manner and going in increasing order of their integer
values. In the first variant, after a forced variable is assigned, the next
variable to branch on is simply the earliest unassigned input in the forward
ordering.

  * $F_{2}$

The default ordering is still “forward,” but now after a forced value is
assigned, the next variable to branch on is the next one _after_ the forced
variable in the forward ordering.

  * $B_{1}$

Similar to strategy $F_{1}$, except the default ordering is “backwards”,
starting with $P(11\dotsb 1)$ and ending with $P(00\dotsb 0)$.

  * $B_{2}$

Similar to strategy $F_{2}$, but with the default ordering being backwards.

For example, in strategy $F_{1}$, the algorithm begins by branching on
$P(000)$. If $P(100)$ were then forced (and there were no variables forced
immediately after), then the next projection to branch on would be $P(001)$,
since $001$ is the first unassigned input in the forward ordering. In
contrast, in the same scenario in $F_{2}$, the algorithm would branch next on
$P(101)$, since $101$ is the next input after the forced input $100$.

###  III.5 Rules 0 and 255

Experimentally, we noticed that coarse-graining to rules 0 or 255 took a very
long time, because pruning could not occur until very late in the search tree.
Note that a nontrivial coarse-graining to (say) rule 0 is the same as a
projection $P$ that is $0$ everywhere except possibly on some input $x$ that
_never occurs_ in the run of the CA $A$, that is, an input $x\in S_{A}^{N}$
that is not in the image of $f_{A}^{N}$. So for coarse-graining to rules 0 or
255, we instead search directly for such non-occuring inputs $x\in S_{A}^{N}$.
When such an input is found, if coarse-graining to rule 0 we simply set
$P(x)=1$, to ensure our projection is not the trivial all-0s projection, which
recall we are excluding by fiat (for rule 255 we set $P(x)=0$ to avoid the
all-1s projection).

If every possible projection input occurs at least once in the image of
$f_{A}^{N}$, then we immediately conclude that $A$ cannot be non-trivially
coarse-grained to either rule 0 or 255. Also, note that by this
characterization, a rule $A$ can be non-trivially coarse-grained to rule 0 if
and only if it can be coarse-grained to rule 255. (Note that for most $A$,
this does _not_ follow from the fact that $\sigma_{01}(0)=255$, since that
would only let us conclude that $A\to 0\Leftrightarrow\sigma_{01}(A)\to 255$.)

##  IV Results

###  IV.1 Coarse grainings

By implementing our improved method, we significantly sped up the search for
coarse-graining. We used our new method(s) to exhaustively find all coarse-
grainings between elementary CA with supercell size $N$ up to 7, on a
commodity laptop. Israeli & Goldenfeld reported results up to $N=4$, and
briefly discussed further results they achieved with large amounts of time on
a super-computer. Our Figure 5 extends their results up to $N=7$.

Figure 5:  Shows the coarse-graining transitions within the 256 elementary CA.
Results from supercell size $N=2$ to $N=7$. An arrow indicates that the first
rule may be coarse grained to the second rule for at least one choice of
supercell size $N$ and nontrivial projection $P$.

The following coarse-grainings were discovered only at supercell sizes 5 or 6,
but not smaller (CA are bracketed according to their symmetry classes):

  * •

The following rules coarse-grain to rules 0 and 255: [25, 61, 67, 103], [41,
97, 107, 121], [43, 113], [54, 147], [57, 99], [62, 118, 131, 145], [73, 109],
[94, 133], [104, 233], [110, 124, 137, 193], [122, 161], [142, 212],

  * •

The following coarse-grain to rule 204: 23, [36, 219], [50, 179], 77, [108,
201], [132, 222], 178, 232

  * •

$[1,127],[19,55]\to 51$

  * •

$[22,151],[104,233]\to[128,254]$

  * •

$[66,189]\to 170$, and the symmetric $[24,231]\to 240$

After $N=5,6$, despite exhaustively and conclusively searching at $N=7$, we
did not discover any new pairs $(A,B)$ such that $A$ coarse-grained to $B$ at
supercell size $7$ but not smaller.

The main diagram of results in [6, 7] has 54 elementary CA that have no
nontrivial coarse-grainings at all up to supercell size 4; using a super-
computer to explore larger supercell sizes, Israeli & Goldenfeld then report
that there were only 16 elementary CA for which they couldn’t find a
nontrivial coarse-graining. On a laptop, we find only 20 elementary CA that do
not have coarse-grainings up to supercell size $N=7$: the 16 previously
reported, along with the two symmetry classes [37,91] and [164, 218]. The non-
trivial coarse-grainings for these four are presumably present at larger
supercell sizes.

###  IV.2 Comparison of algorithms

We compared the four implementation variants $F_{1},F_{2},B_{1},B_{2}$ of our
method discussed above with the brute force method of [6, 7]. The brute force
approach does not include any of the new ideas developed in this paper, such
as the search tree or pruning branches. The brute force method does, however,
include the reduction due to symmetries.

To measure efficiency we will display the number of comparisons between
$res_{1}$ and $res_{2}$; these are displayed in Figure 6; since the four
graphs for $F_{1},F_{2},B_{1},B_{2}$ were so similar, we also report the
numerical values for the five approaches in Table 1.

Figure 6:  Number of comparisons $res_{1}=res_{2}$ used by each of the five methods—brute force, $F_{1},F_{2},B_{1},B_{2}$—in finding coarse-grainings up to supercell size $N=6$. (The plots for $F_{1},B_{1}$ aren’t visible simply because they’re so well-overlapped by those for $F_{2},B_{2}$ at this scale.) Note the $y$-axis is scaled logarithmically, so all four of our new approaches appear to exhibit singly exponential scaling (approximately $\propto 10^{N}$), in comparison to the $2^{2^{N}}$ scaling of the brute force algorithm. The counts for the brute force algorithm are not reported for $N>4$ as the algorithm took too long to run.  $N$ | Brute Force | $F_{1}$ | $F_{2}$ | $B_{1}$ | $B_{2}$
---|---|---|---|---|---
2 | 1,013K | 755K | 755K | 731K | 731K
3 | 28,942K | 8,169K | 8,165K | 7,990K | 8002K
4 | 16,687M | 78,964K | 78,830K | 79,954K | 77,999K
5 | * | 694M | 694M | 685M | 686M
6 | * | 6,778M | 6,722M | 5,777M | 5,798M
Table 1:  The number of comparisons of $res_{1}$ and $res_{2}$ for five
methods (brute force, $F_{1}$, $F_{2}$, $B_{1}$, $B_{2}$). The * indicates
values left uncomputed because the computations took too long.

As expected, the brute force approach was significantly less efficient than
the other 4 selection methods. In Figure 6 all four of our new approaches
appear to exhibit singly exponential $\exp(cN)$ scaling, in comparison to the
$2^{2^{N}}$ scaling of the brute force algorithm. As we progress, pruning
branches apparently has greater effect in reducing the number of comparisons.

Among $F_{1}$, $F_{2}$, $B_{1}$, and $B_{2}$, we see that despite $B_{1}$
being the best approach overall up to $N=6$ there is only about a 1-2%
difference in the number of comparisons for $N=4$ and 5. When we move on to
$N=6$, however, the difference increases to about 15%. On examination, we
found that this greater difference was caused by a single case, the coarse
graining from rule 162 to rule 170. If we eliminate this case, $N=6$ would
also have a miniscule 1-2% difference between our four variants, similar to
that observed at smaller values of $N$.

The coarse-graining from rule 162 to rule 170 at $N=6$ is an interesting
example because of the immense difference between computations in the forward
and backward direction. When applying $B_{1}$, there were precisely 500K
comparisons. Meanwhile, for $F_{1}$, there were 916M comparisons, a factor of
approximately 1800. In this specific example, the backward direction yields a
drastic improvement over the forward direction, but this does not hold true
for all other cases.

In the coarse-grainings from rule 154 to rule 170 with $N=4$, we discovered
that it was the opposite situation in which the forward direction was far more
efficient than the backward direction. When using the backward approach,
66,040 total comparisons were needed to complete. Meanwhile, when we tried the
forward approach, there were only 7,128 comparisons, a factor of just under
10.

We also find examples where there is a significant difference between
selection methods 1 and 2. In coarse-graining rule 29 to rule 51 with $N=7$,
method $B_{1}$ used 1,353K comparisons, while $B_{2}$ used only 241K, a nearly
six-fold difference.

More importantly, we discovered that by altering the initial input projection,
we could reduce the number of comparisons by an extremely large amount.
Although it was a coincidence, we found that when choosing to first branch on
$P(1010101)$, only 60 comparisons were required to complete (just 60, not
60K!). Even if we do not expect such drastic improvements in all cases, the
magnitude of this drop was surprising, and suggests that heuristics for better
selection of projection inputs to branch on could yield quite significant
further improvements.

As is to be expected, no one ordering is always the best. Thus, to provide
optimal efficiency we could combine all four of these. Each of the coarse-
grainings are more effective with different types of dynamic orderings due to
their unique binary patterns, so to further reduce the computations, we still
need to discover how to choose a selection method that will be most compatible
with each of the rules. We now discuss a few patterns we noticed that might be
useful in this endeavor.

Using these new approaches, we observed four common patterns in the cuts. The
first in which both the 0 and 1 branches immediately prune off at level 1, the
lowest level. The second pattern tapers from higher to lower levels in the 0
branch but cuts immediately at level 1 in the 1 branch. The third pattern is
the opposite, in which 0 immediately cuts at level 1 but the 1 branch tapers
from lower to higher levels. The last pattern we saw frequently was a
combination of the second and third, in which the 0 branch tapers from higher
to lower levels while the 1 branch tapers from lower to higher levels.

While most coarse grainings followed these four general trends with slight
deviations, there were some outliers. The coarse-graining from rule 162 to
rule 170 in the forward ordering was unique among all those we explored. It
had a seemingly highly irregular pattern of cuts, with most cuts occurring
only very late in the tree, having little effect in reducing the runtime.
However, when we tried reversing the order, starting at 11, the shape
transformed into the fourth pattern mentioned above and the computation
finished fairly quickly.

Overall, we saw that when there were cases in which the coarse-graining did
not follow one of the preceding four patterns, they far took longer to compute
than others. However, as seen in the example above, changing the direction can
sometimes help fix this issue.

##  V Conclusion

We developed a much more efficient algorithm for finding coarse-grainings
between cellular automata, using a backtracking search, with propagation of
forced values and a dynamic variable ordering. Experimentally, our method
appears to have singly-exponential time scaling, compared to the previous
brute-force method whose runtime scaled as $2^{2^{N}}$ [6, 7]. Using our
method we could examine exhaustively all possible coarse-grainings of
elementary cellular automata up to supercell size 7, extending the previous
results up that only went up to supercell size 4 [6, 7]. We found 26 new
symmetry classes of coarse-grainings (56 coarse grainings total).

To explore further improvements to our algorithm, we examined several
different variable orderings, and found that each had different advantages. We
examined specific cases that suggest that heuristics to select which variable
to branch on next could have quite drastic effects on the runtime. Other
possible improvements are suggested by analogy with the SAT literature, such
as an analogue of clause learning, or to try to incorporate SAT solvers
directly but dynamically (rather than an all-at-once reduction, see Remark 1).

Interesting open questions include the original questions raised in [6,
7]—such as whether the remaining 16 elementary CA can be coarse-grained at
all—as well as new questions we highlight, such as classifying the complexity
of the coarse-graining problem. The experimentally observed singly-exponential
scaling of our algorithm raises the possibility that the problem might in fact
be in the complexity class $\mathsf{EXP}$, making the question of whether it
is in $\mathsf{EXP}$ or is $\mathsf{NEXP}$-complete even more salient. There
also remains the question finding mathematical methods to prove when there
does not exist a coarse graining between two given (elementary) cellular
automata— _regardless_ of supercell size—which would enable us to complete the
picture of coarse-grainings between elementary CA.

###### Acknowledgements.

We would like to thank ATHENA by WiSTEM, the organization that matched the
authors together, eventually leading to this project. The authors were
partially funded by NSF grant DMS-1829826 (formerly DMS-1622390).

## References

  * von Neumann [1966] J. von Neumann,  _Theory of self-reproducing automata_ (Univeristy of Illinois Press, 1966) edited and completed by A. W. Burks.
  * B. and Edelstein-Keshet [1993] E. G. B. and L. Edelstein-Keshet, Cellular automata approaches to biological modeling, J. Theor. Biol. 160, 97 (1993).
  * Ilachinski [2001] A. Ilachinski,  _Cellular automata_ (World Scientific Publishing Co., Inc., River Edge, NJ, 2001) a discrete universe.
  * Mitchell [1998] M. Mitchell, Computation in cellular automata: A selected review, in  _Nonstandard Computation_ (Weinheim: VCH Verlagsgesellschaft., 1998) pp. 95–140.
  * Sarkar [2000] P. Sarkar, A brief history of cellular automata, ACM Comput. Surv. 32, 80 (2000).
  * Israeli and Goldenfeld [2004] N. Israeli and N. Goldenfeld, Computational irreducibility and the predictability of complex physical systems, Phys. Rev. Lett. 92, 074105 (2004).
  * Israeli and Goldenfeld [2006] N. Israeli and N. Goldenfeld, Coarse-graining of cellular automata, emergence, and the predictability of complex systems, Phys. Rev. E 73, 026203 (2006).
  * Wolfram [2002] S. Wolfram,  _A new kind of science_ (Wolfram Media, Inc., Champaign, IL, 2002).
  * Wolfram [1984a] S. Wolfram, Universality and complexity in cellular automata, Phys. D 10, 1 (1984a), cellular automata (Los Alamos, N.M., 1983).
  * Wolfram [1984b] S. Wolfram, Cellular automata as models of complexity, Nature 311, 419 (1984b).
  * Davis and Putnam [1960] M. Davis and H. Putnam, A computing procedure for quantification theory, J. Assoc. Comput. Mach. 7, 201 (1960).
  * Davis  _et al._ [1962] M. Davis, G. Logemann, and D. Loveland, A machine program for theorem-proving, Commun. Assoc. Comput. Mach. 5, 394 (1962).
  * Magiera and Dzwinel [2014] K. Magiera and W. Dzwinel, A novel algorithm for coarse-graining of cellular automata, in  _Cellular Automata - 11th International Conference on Cellular Automata for Research and Industry, ACRI 2014, Krakow, Poland, September 22-25, 2014. Proceedings_, Lecture Notes in Computer Science, Vol. 8751, edited by J. Was, G. C. Sirakoulis, and S. Bandini (Springer, 2014) pp. 258–267.
  * Dzwinel and Magiera [2015] W. Dzwinel and K. Magiera, Irreducible elementary cellular automata found, J. Comput. Sci. , 300 (2015).


# Theory and Computation of Substructure Characteristic Modes

Mats Gustafsson, , Lukas Jelinek, Miloslav Capek, ,
Johan Lundgren, , and Kurt Schab Manuscript received ; revised . This work was
supported by ELLIIT - an Excellence Center at Linkoping-Lund in Information
Technology and the Czech Science Foundation under project No. 21-19025M.M.
Gustafsson and J. Lundgren are with Lund University, Lund, Sweden, (e-mails:
{mats.gustafsson, johan.lundgren}@eit.lth.se).L. Jelinek and M. Capek are with
Czech Technical University in Prague, Czech Republic (e-mails:
{lukas.jelinek,miloslav.capek}@fel.cvut.cz).K. Schab is with the Santa Clara
University, Santa Clara, USA (e-mail: kschab@scu.edu).Color versions of one or
more of the figures in this paper are available online at
http://ieeexplore.ieee.org.

###### Abstract

The problem of substructure characteristic modes is reformulated using a
scattering matrix-based formulation, generalizing subregion characteristic
mode decomposition to arbitrary computational tools. It is shown that the
scattering formulation is identical to the classical formulation based on the
background Green’s function for lossless systems. The scattering formulation,
however, opens a variety of new subregion scenarios unavailable within
previous formulations, including cases with lumped or wave ports or subregions
in circuits. Thanks to its scattering nature, the formulation is solver-
agnostic with the possibility to utilize an arbitrary full-wave method.

###### Index Terms:

Antenna theory, characteristic modes, computational electromagnetics,
eigenvalues and eigenfunctions, scattering.

##  I Introduction

Characteristic mode decomposition [1, 2, 3] plays an important role in the
design [4] of antennas, such as electrically small antennas, MIMO systems, and
arrays. A recent extension of the scattering-based formulation of
characteristic mode decomposition [2] has broadened its application scope to
include arbitrary electromagnetic solvers [5, 6, 7], enabling advanced
applications with arbitrary material distributions. Despite its numerous
advantages [7] and rapid evaluation capabilities [8], the scattering approach
to characteristic modes was not readily extended to the substructure variant
[9] frequently applied using impedance-based methods [3].

Substructure characteristic mode decomposition involves dividing the
scattering scenario into a controllable (or accessible, see [9, 10]) region
and a background. In this approach, any structural modification, such as
antenna design or selective excitation, is confined to the controllable
region. An example is a patch antenna situated over a ground plane, where the
ground plane has a significant impact on most characteristic modes, yet the
designer can effectively influence only the patch design region [11, 9, 12].
The substructure characteristic mode decomposition focuses on modes most
closely associated with the controllable region by using altered forms of
operators describing the scattering problem. This makes substructure modes an
attractive approach for studying the behavior of radiating devices affected by
nearby objects like vehicles, electronic platforms, or biological tissues
[13].

Following the state-of-the-art procedures, substructure characteristic modes
can be interpreted as a classical formulation of characteristic mode
decomposition where the underlying Green’s dyadic includes the surroundings
(background) of the studied region [14, 15, 16, 17]. One implementation of
this approach using the method of moments (MoM) [18] formulations is based on
a Schur complement [9, 19], which constructs a compressed impedance matrix for
a scatterer in the presence of background objects.

This paper expands the definition of substructure modes by introducing a
scattering-based variant, formulated through a generalized eigenvalue problem
of two scattering matrices [20, § 7.8.1], [21, § 4.3]. The first matrix
accounts for the complete scattering problem, while the second represents the
background. Following assumptions widely used in characteristic mode analysis,
we assume losslessness, linearity, and time-invariance for all materials. The
employed scattering matrices can be easily substituted with transition
matrices (T-matrices) [22] or by scattering dyadics [20, § 4.3], which can be
constructed by an arbitrary numerical technique [5, 6, 7]. The iterative
algorithm [8] adapted here for the evaluation of substructure modes mitigates
the computational burden associated with increasing electrical size and model
complexity. This makes the proposed formulation applicable to various problems
typically outside the scope of characteristic mode analysis, including
microwave circuits and optical circuits spanning large ranges of electrical
sizes with complex material distributions.

The paper is organized as follows. Section II briefly recapitulates the
meaning of substructure characteristic modes and the used nomenclature.
Section III introduces a novel scattering-based formulation. Its equivalence
to MoM-based substructure characteristic mode decomposition is shown in
Section IV on cases when the formulation was known before, including perfect
electric conductors, a hybrid of the method of moments and T-matrix method,
and infinite ground plane. The advantages of scattering formulation are shown
in Section V, dealing with challenging examples that are impractical to
compute using previous approaches. The outcomes are discussed in Section VI,
and the paper is concluded in Section VII.

##  II Substructure Characteristic Modes

Insights into essential features of substructure characteristic modes can be
obtained by comparing Fig. 1 and Fig. 2, where we adapt a PIFA-like example
[9] made of perfect electric conductor (PEC), illustrating the characteristic
modal significance of the entire region
$\varOmega=\varOmega_{\mathrm{c}}\cup\varOmega_{\mathrm{b}}$ and substructure
modal significance corresponding to the controllable region
$\varOmega_{\mathrm{c}}$ in the presence of the background region
$\varOmega_{\mathrm{b}}$.

In Fig. 1, the entire device is studied, and the resulting modal spectrum
includes modes which can be induced and superimposed by characteristic
excitations applied over the entire structure. This results in a dense set of
modes with high modal significance. However, when excitations are confined to
the design region $\varOmega_{\mathrm{c}}$, treating the region
$\varOmega_{\mathrm{b}}$ as a fixed background, the resulting modal spectrum,
shown in Fig. 2, exhibits notable differences. Here, only modes that can be
induced or superimposed using excitations on the design region
$\varOmega_{\mathrm{c}}$ are present, resulting in a much sparser
eigenspectrum and revealing a resonance peak around $1.8\,\mathrm{GHz}$ which
is supported by currents on the PIFA. It is important to note that, in this
case, modal current distributions exist over both regions
$\varOmega_{\mathrm{c}}$ and $\varOmega_{\mathrm{b}}$. Nevertheless, the
absence of independent control of excitations over the region
$\varOmega_{\mathrm{b}}$ effectively removes currents over that region from
the set of degrees of freedom defining each mode.

Figure 1: Modal significances $|t_{n}|$ computed for a PEC geometry adapted
from [9]. The dimensions read $\ell=120\,\mathrm{mm}$, $w=60\,\mathrm{mm}$,
$h=15\,\mathrm{mm}$, and $d=30\,\mathrm{mm}$. Impedance-based formulation [3]
and scattering-based formulation (1) are employed with
$\mathbf{S}_{\mathrm{b}}$ being the unit matrix. Figure 2: Modal significances
$|t_{n}|$ computed for the substructure $\varOmega_{\mathrm{c}}$ in the
presence of the background scatterer $\varOmega_{\mathrm{b}}$. Dimensions are
identical to Fig. 1. The impedance-based method relies on formulation [9],
while the scattering-based formulation relies on (1).

Previously, calculations like these were limited by their reliance on MoM [9,
23]. This paper demonstrates that substructure characteristic modes can be
evaluated using the scattering matrix formalism, significantly expanding their
scope and computational capabilities.

##  III Scattering Formulation of Substructure Characteristic Modes

Assume a controllable antenna region $\varOmega_{\mathrm{c}}$ and the
surrounding (or background) region $\varOmega_{\mathrm{b}}$, as illustrated in
Fig. 3. Let $\mathbf{S}_{\mathrm{b}}$ denote that scattering matrix [21, 20]
for the background region $\varOmega_{\mathrm{b}}$, which connects incoming
waves represented by coefficients collected in a vector $\mathbf{a}$ into
outgoing waves represented by a vector $\mathbf{f}$. Similarly, let
$\mathbf{S}$ denote the scattering matrix for the composite object
$\varOmega_{\mathrm{c}}\bigcup\varOmega_{\mathrm{b}}$. The core hypothesis of
this paper is that substructure characteristic modes are determined from the
generalized eigenvalue problem

$\mathbf{S}\mathbf{a}_{n}=s_{n}\mathbf{S}_{\mathrm{b}}\mathbf{a}_{n},$ (1)

where $\mathbf{a}_{n}$ are characteristic excitations and $s_{n}$
characteristic scattering eigenvalues. Throughout the remainder of the paper,
we refer to these quantities as characteristic modes due to their equivalence
with the characteristic currents $\mathbf{I}_{n}$ typically used in impedance-
based formulations [3, 5]. This formulation generalizes the one defined by
[24] and developed in [5] for objects in free space, which implicitly
considers the background scattering matrix $\mathbf{S}_{\mathrm{b}}$ to be an
identity matrix. It also covers the special case studied for periodic
structures in [25], where the background scattering matrix
$\mathbf{S}_{\mathrm{b}}$ is related to an implicit connection between plane
waves on either side of a scattering surface, analogous to the S-parameters of
a transmission line network. The scattering eigenvalues $s_{n}$ are related to
modal significance $|t_{n}|$ and characteristic eigenvalues $\lambda_{n}$ as
[5]

$t_{n}=\frac{s_{n}-1}{2}\quad\text{and
}\lambda_{n}=-\operatorname{Im}\\{t_{n}^{-1}\\}=\mathrm{j}\frac{s_{n}+1}{s_{n}-1}.$
(2)

Similarly to characteristic modes of isolated objects [5], for lossless
objects, characteristic modes exhibit the following orthogonality relations

$\mathbf{a}_{m}^{\mathrm{H}}\mathbf{a}_{n}=\delta_{mn},\quad\mathbf{f}_{m}^{\mathrm{H}}\mathbf{f}_{n}=\delta_{mn},$
(3)

which, for example, translates to the orthogonality of characteristic far
fields and where H denotes Hermitian transpose and $\delta_{mn}$ is the
Kronecker delta. Equivalence between scattering and impedance formulations is
supported by data shown in Figs. 1 and 2 comparing results from modal
significances produced by (1) and an impedance-based formulation of
substructure modes [9]. A mathematical derivation of this equivalence is
outlined in Section IV and Appendix B.

a)$\varOmega_{\mathrm{b}}$$\mathbf{a}$$\mathbf{f}=\mathbf{S}_{\mathrm{b}}\mathbf{a}$b)$\varOmega_{\mathrm{c}}$$\varOmega_{\mathrm{b}}$$\mathbf{a}$$\mathbf{f}=\mathbf{S}\mathbf{a}$
Figure 3: Sketch of the physical interpretation of scattering matrices
$\mathbf{S}_{\mathrm{b}}$ and $\mathbf{S}$.

###  III-A Variations

For lossless objects, the eigenvalue problem in (1) can be rearranged into
several alternative forms. Scattering matrices are unitary for lossless
objects [21], $\mathbf{S}^{\mathrm{H}}\mathbf{S}=\mathbf{1}$, leading to

$\mathbf{S}_{\mathrm{b}}^{\mathrm{H}}\mathbf{S}\mathbf{a}_{n}=s_{n}\mathbf{a}_{n}\quad\text{and}\quad\mathbf{S}\mathbf{S}^{\mathrm{H}}_{\mathrm{b}}\mathbf{f}_{n}=s_{n}\mathbf{f}_{n},$
(4)

where
$\mathbf{f}_{n}=\mathbf{S}_{\mathrm{b}}\mathbf{a}_{n}=s_{n}^{\ast}\mathbf{S}\mathbf{a}_{n}$
represents the scattered field, $\mathbf{1}$ is the identity matrix, and ∗
denotes complex conjugate. The two versions in (4) are equivalent and differ
only in expressing the eigenvalue problem in incoming (excitation),
$\mathbf{a}_{n}$, or outgoing (scattered or radiated), $\mathbf{f}_{n}$,
waves. We note that the incident field $\mathbf{a}_{n}$ differ from the
scattered field $\mathbf{f}_{n}$ except for the free-space case with
$\mathbf{S}_{\mathrm{b}}=\mathbf{1}$ [5].

Formulation (4) can be implemented using transition matrices111Alternatively,
scattering dyadics [20] can be used in place of the transition matrices
$\mathbf{T}$ and $\mathbf{T}_{\mathrm{b}}$ with no further changes to the
formulation presented in the remainder of this section, save for the
understanding that the transition matrix maps regular spherical waves to
outgoing spherical waves, while the scattering dyadic maps incoming plane
waves to the scattered far field. See [7] for details on these two operators
and their use in characteristic modes. [20]

$\mathbf{S}=2\mathbf{T}+\mathbf{1}\quad\text{and}\quad\mathbf{S}_{\mathrm{b}}=2\mathbf{T}_{\mathrm{b}}+\mathbf{1},$
(5)

which allows for the analysis of substructure characteristic modes solely in
terms of transition matrices and characteristic number $t_{n}$. These
formulations read

$(2\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}\mathbf{T}+\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T})\mathbf{a}_{n}=t_{n}\mathbf{a}_{n}$
(6)

and

$(2\mathbf{T}\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T})\mathbf{f}_{n}=t_{n}\mathbf{f}_{n}$
(7)

and does not follow a scheme analogous to the characteristic modes of isolated
objects [5]. In contrast to (1), relations (4), (6), and (7) enable iterative
matrix-free evaluation [8], an important factor when employing generic
electromagnetic solvers for evaluating substructure characteristic modes for
large problems. Details of the iterative solution is found in Appendix A.
Moreover, relation (6) allows for interpreting the substructure scattered
power as a power of the difference between the scattered field of the
composite object $\mathbf{T}\mathbf{a}_{n}$ and the scattered field of the
background $\mathbf{T}_{\mathrm{b}}\mathbf{a}_{n}$, i.e.

$\frac{1}{2}|(\mathbf{T}-\mathbf{T}_{\mathrm{b}})\mathbf{a}_{n}|^{2}=\frac{-\operatorname{Re}\\{t_{n}\\}}{2}|\mathbf{a}_{n}|^{2}=\frac{1}{2}|t_{n}|^{2}|\mathbf{a}_{n}|^{2}$
(8)

with the interpretation of no substructure scattering for $t_{n}=0$ (or
$s_{n}=1$) and maximum scattering for $t_{n}=s_{n}=-1$ similar to the free-
space case [5]. Here, we have utilized the property [5]
$\mathbf{T}^{\mathrm{H}}\mathbf{T}=-\mathrm{Re}\\{\mathbf{T}\\}$.

##  IV Equivalence Between MoM and Scattering-Based Substructure
Characteristic modes

The proposed formula (1) is equivalent to the impedance-based substructure
modes [9], wherein a structure is separated into two (possibly connected)
regions $\varOmega_{\mathrm{b}}$ (“background”) and $\varOmega_{\mathrm{c}}$
(“controllable”) and excitation is restricted to only the controllable region
$\varOmega_{\mathrm{c}}$. A general scattering problem in the standard
notation for electric field integral equation (EFIE)-formulations [18] reads

$\mathbf{V}=\mathbf{Z}\mathbf{I},$ (9)

where $\mathbf{V}$ is an excitation vector, $\mathbf{Z}$ is the impedance
matrix associated with the scatterer $\varOmega$, and $\mathbf{I}$ is an
induced current density, all represented in a particular basis. Bifurcating
the system into subsets of the basis functions associated with the two regions
$\varOmega_{\mathrm{b}}$ and $\varOmega_{\mathrm{c}}$ and enforcing zero
excitation on the background region leads to [14]

$\mathbf{Z}\mathbf{I}=\begin{bmatrix}\mathbf{Z}_{\mathrm{b}\mathrm{b}}&\mathbf{Z}_{\mathrm{b}\mathrm{c}}\\\
\mathbf{Z}_{\mathrm{c}\mathrm{b}}&\mathbf{Z}_{\mathrm{c}\mathrm{c}}\end{bmatrix}\begin{bmatrix}\mathbf{I}_{\mathrm{b}}\\\
\mathbf{I}_{\mathrm{c}}\end{bmatrix}=\begin{bmatrix}\mathbf{0}\\\
\mathbf{V}_{\mathrm{c}}\end{bmatrix}.$ (10)

By reducing the MoM system $\mathbf{Z}=\mathbf{R}+\mathrm{j}\mathbf{X}$ to its
Schur complement, the substructure characteristic modes are defined as [9]

$\widetilde{\mathbf{X}}\mathbf{I}_{\mathrm{c}n}=\lambda_{n}\widetilde{\mathbf{R}}\mathbf{I}_{\mathrm{c}n}\quad\text{or}\quad\widetilde{\mathbf{Z}}\mathbf{I}_{\mathrm{c}n}=(1+\mathrm{j}\lambda_{n})\widetilde{\mathbf{R}}\mathbf{I}_{\mathrm{c}n},$
(11)

where
$\widetilde{\mathbf{Z}}=\widetilde{\mathbf{R}}+\mathrm{j}\widetilde{\mathbf{X}}=\mathbf{Z}_{\mathrm{c}\mathrm{c}}-\mathbf{Z}_{\mathrm{c}\mathrm{b}}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}$
is an impedance matrix related to the problem-specific numerical Green’s
function [14]. This numerical Green’s function captures the field-generating
behavior of currents on the controllable region $\varOmega_{\mathrm{c}}$ in
the presence of the “background” region $\varOmega_{\mathrm{b}}$.

The proof of equality of (11) and (1) starts with factorization of the
radiation matrix [26] $\mathbf{R}=\mathbf{U}_{1}^{\mathrm{T}}\mathbf{U}_{1}$,
where the matrix $\mathbf{U}_{1}$ projects MoM basis functions onto spherical
waves. Next step is partitioning
$\mathbf{U}_{1}=\begin{bmatrix}\mathbf{U}_{1\mathrm{b}}&\mathbf{U}_{1\mathrm{c}}\end{bmatrix}$
and factorization

$\widetilde{\mathbf{R}}=\widetilde{\mathbf{U}}_{1}^{\mathrm{H}}\widetilde{\mathbf{U}}_{1}$
(12)

for the substructure case, where
$\widetilde{\mathbf{U}}_{1}=\mathbf{U}_{1\mathrm{c}}-\mathbf{U}_{1\mathrm{b}}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}$.
A substructure eigenvalue problem for T-matrix

$\widetilde{\mathbf{T}}\mathbf{f}_{n}=t_{n}\mathbf{f}_{n}$ (13)

with $t_{n}=\left(1+\mathrm{j}\lambda_{n}\right)^{-1}$ can be formulated by
substituting (12) into (11), left multiplication with
$\widetilde{\mathbf{U}}_{1}\widetilde{\mathbf{Z}}^{-1}$ and by identifying

$\widetilde{\mathbf{T}}=-\widetilde{\mathbf{U}}_{1}\widetilde{\mathbf{Z}}^{-1}\widetilde{\mathbf{U}}_{1}^{\mathrm{H}}\quad\mathrm{and}\quad\mathbf{f}_{n}=-\widetilde{\mathbf{U}}_{1}\mathbf{I}_{\mathrm{c}n}.$
(14)

Appendix B then shows that matrix $\widetilde{\mathbf{T}}$ equals to the
matrix used in (7). Similarly to full characteristic modes, the characteristic
substructure current can be evaluated as [5]

$\mathbf{I}_{\mathrm{c}n}=t_{n}^{-1}\widetilde{\mathbf{Z}}^{-1}\widetilde{\mathbf{U}}_{1}^{\mathrm{H}}\mathbf{f}_{n}.$
(15)

###  IV-A Example — Demonstrating equivalence

To demonstrate the numerical equivalence between the impedance and scattering
formulations in (11) and (1), we employ a PEC structure reported in [9], which
consists of a PIFA-like region connected to a finite ground plane, as depicted
in the inset of Fig. 1. In Fig. 1, the characteristic modes of the full
structure are computed without any separation into controllable and background
regions, i.e., the entire region is controllable. Impedance and scattering
operators are calculated using MoM, and the modal significances produced by
the impedance and scattering formulations are numerically identical.
Similarly, Fig. 2 displays modal significances when only part
$\varOmega_{\mathrm{c}}$ of the structure is controllable (following [9]), and
again the impedance (11) and scattering formulations (1) agree within
numerical precision.

###  IV-B MoM & T-matrix hybrid

Hybridization of MoM and T-matrix techniques efficiently models complex
inhomogeneous structures, so long as the regions described by each method can
be separated by a plane [27]. Transition matrices can be computed for
arbitrary background problems, but the computational acceleration over full
MoM implementations is most pronounced when analytic forms of transition
matrices can be employed, e.g., Mie series results for layered spherical
structures or infinite ground planes. In this section, we discuss how this
form of hybridization also allows for the efficient computation of
substructure modes when part or all of the background region is represented by
a transition matrix.

Consider a case from Fig. 2 to which a dielectric object is introduced as in
Fig. 4. The scattering properties of a dielectric object can be efficiently
described by matrix $\mathbf{T}_{\mathrm{\mathrm{b}1}}$, while the scattering
properties of the metallic structure are described by the impedance matrix
$\mathbf{Z}$ composed according to (10).

Figure 4: Modal significances $|t_{n}|$ computed for the substructure
$\varOmega_{\mathrm{c}}$ in the presence of the background which contains
scatterer $\varOmega_{\mathrm{b}}$ and a dielectric sphere filled with
relative permittivity $\varepsilon_{\mathrm{r}}=4$. The diameter of the sphere
equals 180 mm and the distance between the sphere and the metallic structure
is 30 mm. Other dimensions are identical to Fig. 1. The impedance-based method
relies on (11) and (16), while the scattering-based formulation relies on (1).

The reaction of the entire setup to an external excitation
$\mathbf{a}^{\mathrm{i}}$ can be written as [27]

$(\mathbf{Z}+\mathbf{U}^{\mathrm{T}}_{4}\mathbf{T}_{\mathrm{\mathrm{b}1}}\mathbf{U}_{4})\mathbf{I}=(\mathbf{U}_{1}^{\mathrm{T}}+\mathbf{U}_{4}^{\mathrm{T}}\mathbf{T}_{\mathrm{\mathrm{b}1}})\mathbf{a}^{\mathrm{i}}.$
(16)

The matrix
$\mathbf{Z}+\mathbf{U}^{\mathrm{T}}_{4}\mathbf{T}_{\mathrm{\mathrm{b}1}}\mathbf{U}_{4}$
is the impedance matrix of a metallic body in the presence of a dielectric
object, with the second term being interpreted as a contribution of matrix
$\mathbf{T}_{\mathrm{\mathrm{b}1}}$ to the background Green’s function [16].
Characteristic decomposition of this modified impedance matrix gives a
spectrum similar to Fig. 1, with the entire conducting region in the presence
of a dielectric sphere [28]. On the other hand, if substructure modes solely
excitable from region $\varOmega_{\mathrm{c}}$ are desired, then this matrix
is further modified and decomposed according to (11). The resulting spectrum
of modal significances, where the entire region $\varOmega_{\mathrm{b}}$ is
considered a background, is shown in Fig. 4. The equivalence of this
formulation to (1) is demonstrated in the same figure.

###  IV-C Infinite ground plane

Figure 5: Modal significances $|t_{n}|$ computed for the substructure
$\varOmega_{\mathrm{c}}$ in the presence of the background
$\varOmega_{\mathrm{b}}$, which in this case is an infinite perfectly
conducting ground plane. Dimensions are identical to Fig. 1.

An extreme version of the structure studied in Fig. 2 is a region
$\varOmega_{\mathrm{c}}$ placed above an infinitely large perfectly conducting
plane as depicted in Fig. 5. When substructure characteristic modes are
evaluated in this scenario, a common practice is to construct and decompose
the matrix $\widetilde{\mathbf{Z}}_{\mathrm{c}\mathrm{c}}$ for the region
$\varOmega_{\mathrm{c}}$ using the Green’s function based on equivalent image
currents [11]. Another possibility is to decompose the impedance matrix
belonging to region $\varOmega_{\mathrm{c}}$ and its mirror image and then
employ point symmetries [29] to filter out modes belonging to irreducible
representation with even parity (those would belong to a perfectly
magnetically conducting ground plane).

Within the scattering formulation (1), the problem is solved by adding an
image of region $\varOmega_{\mathrm{c}}$ and an image of the incident field.
This results in a total electric field with vanishing tangential components at
the ground plane $\varOmega_{\mathrm{b}}$. The symmetry of the incident field
eliminates half of the spherical waves. Specifically, only spherical waves
exhibiting electric dipole moment normal to the ground plane and tangential
magnetic dipole moment remain. For example, placing the ground plane
$\varOmega_{\mathrm{b}}$ in the $xy$ plane and denoting $l,m$ the degree and
azimuthal numbers, respectively, only TE spherical waves with even $l+m$ and
TM spherical waves with odd $l+m$ remain. The equivalence of the impedance-
based formulation and scattering formulation is presented in Fig. 5.

The practicality of using substructure characteristic modes is demonstrated by
comparing spectra from Figs. 2, 4 and 5. In all cases, the black resonance
peak appearing between $1.5\,\mathrm{GHz}$ and $2\,\mathrm{GHz}$ can be
identified but would be lost if full characteristic decomposition was made,
cf., Fig. 1. This resonance belongs to a mode responsible for the PIFA
operation, and substructure decompositions show how it is affected by a
particular background.

##  V Utilization of Arbitrary Full-Wave Solver

In case no in-house code is available to construct scattering matrices
$\mathbf{S}$ and $\mathbf{S}_{\mathrm{b}}$, commercial simulators can be
employed instead. This requires an interface between the simulator and a post-
processor to assemble and decompose the scattering matrices.

To demonstrate the flexibility of scattering formulation (1), the scattering
dyadic matrices [7] are employed in this section to evaluate sub-structure
characteristic modes on several complex examples using a commercial solver as
the core computational engine. To mitigate the computational burden stemming
from the fact that the full-wave evaluation is repetitively performed for all
columns of the scattering dyadic matrices, an iterative procedure [8] is
employed and modified for the substructure case, see Appendix A. Codes for
these examples implemented in MATLAB and Altair FEKO [30] are available at
[31].

###  V-A Planar Inverted-F Antennas

A PIFA of the same dimensions as in Fig. 1 is studied once more, now with the
addition of two example configurations of the dielectric substrate. Modal
significance data for these examples, along with the previously studied PEC-
only model, are shown in Fig. 6.

The reference case (PIFA #1) is made solely of PEC and represented by black
traces. This case serves only as a verification of the scattering dyadic
matrix procedure [7]. The resulting data are indistinguishable from the curves
in Fig. 2. The two cases involving dielectrics are also shown in Fig. 6 and
include two configurations of a lossless dielectric material with relative
permittivity $\varepsilon_{\mathrm{r}}=2.33$. One case (PIFA #2) is built with
two thin (1.575 mm) dielectric layers backing the parallel conductive layers
in the PIFA system. The other case (PIFA #3) includes a dielectric block
spanning the entire space between the ground plane and the upper motif.

It can be seen that the eigen-traces of PIFA #1 and PIFA #2 are comparable,
which is given by the fact that thin substrate minimizes the effect of the
dielectrics, i.e., the effective permittivity between the ground plane and the
upper motif is close to unity, leading to only a slight downward shift in the
resonant frequency of the dominant characteristic mode. This is not the case
with PIFA #3, which behaves both qualitatively and quantitatively differently.
The modal significance maxima are considerably shifted towards lower
frequencies due to the reduced wavelength within the dielectric. The shape of
the trace associated with the dominant mode is also visually quite different.

Figure 6: The first three substructure characteristic modes for different
planar F antennas. The first model is taken from Fig. 2 and serves as a
verification. The second model is built on a substrate with a thickness of
1.575 mm and relative permittivity $\varepsilon_{\mathrm{r}}=2.33$. The third
model has a dielectric block of $\varepsilon_{\mathrm{r}}=2.33$ filling the
space between the ground plane and the antenna. All models are lossless and
treated with an iterative algorithm that employs surface equivalence method-
of-moment from Altair FEKO [30]. The traces were adaptively refined to obtain
smooth curves with relatively few samples.

###  V-B Simplified CubeSat Model

The previous example demonstrates the quantitative effect of dielectric
substrates on modal characteristics of a PIFA, which can, to a certain extent,
be anticipated from general rules of antenna design, i.e., shifting of
resonances to lower frequencies through the use of dielectric loading. In this
final example, however, we consider the modal characteristics of a complex
system for which limited engineering intuition can be applied a priori.

The substructure characteristic modes of a simplified CubeSat model are
evaluated in Fig. 7. The CubeSat frame of dimensions $10\times 10\times 15$
cm3 (approximately 1.5U format [32]) is formed by a PEC strip of width 1 cm.
The frame is galvanically connected with the solid PEC top cover and left open
on the underside. The interior is filled with a dielectric material with
relative permittivity $\varepsilon_{\mathrm{r}}=3$. The PIFA of dimensions
$w=5$ cm, $d=2.5\,$cm, and $h=1.25\,$cm is mounted on the top of the CubeSat
and, similarly as in Fig. 2, only this antenna region is considered
controllable for modal analysis, see the inset in Fig. 7. Notice the antenna
size was reduced by 20% as compared to Fig. 2 to reflect the change of the
ground plane size (from 120 cm to 100 cm) in this example, following
dimensions of standard CubeSats. The model is lossless and treated with an
iterative algorithm that employs surface equivalence method-of-moment from
Altair FEKO [30]. The traces were adaptively refined to obtain smooth curves
with relatively few samples.

The spectrum shown in Fig. 7 is dominated by only one mode, an observation
agreeing well with the similar arrangement in Fig. 2. Considering the 20%
relative size reduction of PIFA in Fig. 7 as compared with Fig. 2, the
dominant mode resonates at comparable frequency. The background structure
contains a large dielectric block with many internal resonances (approximately
80 cavity resonances between 0.1 GHz and 2.6 GHz). Nevertheless, all these
resonances are filtered out by the substructure formulation (1). The only
exception is the mode having an abrupt increase of modal significance around
1.6 GHz. This spike is still well-modeled by the iterative algorithm by finely
sampling the frequency axis.

Figure 7: The first five substructure characteristic modes for a simplified
model of a CubeSat. The CubeSat has dimensions $10\times 10\times 15$ cm3. The
frame is made of PEC and has a width of 1 cm. The top cover is also made of
PEC and forms a ground plane for a PIFA mounted on the CubeSat. The PIFA has
dimensions of $w=5$ cm, $d=2.5$ cm, $h=1.25$ cm, and is the only controllable
region in the model. The entire CubeSat is filled with dielectrics of relative
permittivity $\varepsilon_{\mathrm{r}}=3$.

##  VI Discussion

The theoretical developments and set of examples provided in previous sections
thoroughly demonstrate that the scattering formulation of substructure
characteristic modes (1) is able to reproduce cases treated using classical
formulation [9], while also allowing for the analysis of arbitrarily complex
material distributions without modifying the evaluation procedure.

With the exception of the CubeSat example shown in Fig. 7, all examples are
based on the same PIFA-like controllable region in the presence of varying
backgrounds. Substantial differences between the resulting modal
characteristics in each example clearly illuminate the high potential impact
of background objects on modal performance. This is particularly noticeable in
Fig. 6, where the inclusion of thin dielectric support layers leads to non-
negligible changes in modal significance. If outputs from characteristic mode
analyses are used in the design of antennas, it is therefore critical to fully
model any background objects, including dielectrics, rather than using
simplified background models and small perturbation approximations. The
approach presented in this manuscript facilitates this rigorous analysis.

To further investigate the potential of the scattering formulation, it is
worth considering the scenarios depicted in Fig. 8. All the cases treated in
Sections II, IV and V solely dealt with panel (a), in which controllable and
background regions might share a boundary but are otherwise disjoint.

Figure 8: Illustrations of different setups to be considered for decomposition
via substructure characteristic modes. (a) Classical arrangement with non-
overlapping (possibly connected) regions. (b) Region $\varOmega_{\mathrm{c}}$
partly or fully overlapping with background region $\varOmega_{\mathrm{b}}$.
(c) Region $\varOmega_{\mathrm{b}}$ in a cavity and $\varOmega_{\mathrm{c}}$
also contains waveguide ports. (d) Region $\varOmega_{\mathrm{c}}$ overlapping
with background material in the presence of ports. Regions
$\varOmega_{\mathrm{c}}$ and $\varOmega_{\mathrm{b}}$ are composed of
arbitrary material distributions.

The first generalization is shown in panel (b), where controllable and
background regions share the same volumes yet have different material
properties. The meaning of such a situation can be understood from the volume
equivalence principle. When building an equivalent description of a given
scattering scenario, one can choose which part of polarization belongs to
equivalent sources and which belongs to a background. In such a case, the
controllable region is not a physical structure but a contrast between two
material distributions. In principle, this scenario can be approached by
classical treatment using a partial background Green’s function, however, the
computation overhead will be considerable. On the other hand, within the
scattering formulation, the problem from panel (b) is treated exactly like the
problem from panel (a), i.e., by separately evaluating scattering from the
entire system and the background.

The second generalization is the addition of ports, a situation depicted in
panel (c). A possible scenario might be a substructure problem of Section II,
where a port is added to controllable degrees of freedom222If the port is
considered part of the background, the system will be lossy, which is not
considered in this paper.. Although proposals exist for evaluating
characteristic modes on antennas loaded by ports [33], we stress here that the
ports present no qualitative change for the scattering formulation.
Furthermore, the scattering formulation suggests that substructure
characteristic modes can be evaluated even for circuits [21]. For example, the
scattering matrices $\mathbf{S}_{\mathrm{b}}$ might represent a circuit
without a particular region, which will be considered the controllable part.

Lastly, panel (d) of Fig. 8 combines all preceding scenarios. A possible
application of this scenario is the analysis and optimization of minimum
scattering antenna [34] by inspecting characteristic modes and their
significance [35], or techniques involving generalized scattering matrix [36,
37]. In such a case, the substructure characteristic modes are not reachable
via the Schur complement method [9], and their formulation will be challenging
if background Green’s function is employed. Yet, the scattering formulation
stays always the same, solely demanding to evaluate scattering matrices
$\mathbf{S}$ and $\mathbf{S}_{\mathrm{b}}$ for a situation with and without
the controllable part of the structure.

##  VII Conclusion

The method proposed in this paper extends the scattering-based formulation of
characteristic modes to general substructure problems. Such problems were
previously limited to cases where numerical or analytical problem-specific
(non-free-space) Green’s functions are known, typically by way of Schur
complement methods based on the method of moments. Using the proposed
extension of scattering-based characteristic mode analysis, this requirement
is lifted, and characteristic modes for all linear problems with arbitrary
structure/substructure designations can be analyzed.

The examples throughout the paper demonstrate the flexibility of the method
through its application to a family of problems based on a PIFA-like antenna
above a finite ground plane. The resulting data demonstrate equivalence
between the proposed method and impedance-based formulations, and highlight
the generality of the approach to variations involving large dielectric
background media, infinite ground planes, and finite dielectric regions.
Further discussion regarding its application to problems involving wave ports
and acceleration using iterative algorithms points to several areas for
continued research, with the general direction aiming toward the fast and
efficient characteristic mode analysis of antennas and scatterers in highly
complex, arbitrary material environments.

##  Appendix A Matrix-Free Evaluation of Substructure Characteristic Modes

Characteristic modes can be evaluated efficiently using an iterative matrix-
free algorithm [8]. Instead of having full, explicit knowledge of a matrix
being decomposed, these algorithms rely on knowing the result of a matrix-
vector multiplication of that matrix with an arbitrary vector. In the case of
characteristic modes, the matrix itself is a scattering operator, while the
matrix-vector multiplication represents the solution to a scattering problem
for a particular excitation. Nevertheless, algorithms such as Arnoldi
iteration used in [8] are not suitable for generalized eigenvalue problems,
such as (1). Therefore, to evaluate substructure characteristic modes in a
matrix-free fashion, the formulations (4), (6) and (7) are used. Furthermore,
it is important to realize that only results of multiplications
$\mathbf{T}\mathbf{x},\mathbf{T}_{\mathrm{b}}\mathbf{x},\mathbf{S}\mathbf{x},\mathbf{S}_{\mathrm{b}}\mathbf{x}$
are accessible via general purpose electromagnetic solvers. When adapting the
algorithms described in [8], the following changes in the desired matrix-
vector products are made

$\displaystyle\mathbf{S}_{\mathrm{b}}^{\mathrm{H}}\mathbf{S}\mathbf{a}_{n}$
$\displaystyle=\left(\mathbf{S}_{\mathrm{b}}\left(\mathbf{S}\mathbf{a}_{n}\right)^{*}\right)^{*}$
(17)
$\displaystyle(2\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}\mathbf{T}+\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T})\mathbf{a}_{n}$
$\displaystyle=\left(\mathbf{T}_{\mathrm{b}}\left(\mathbf{a}_{n}+2\mathbf{T}\mathbf{a}_{n}\right)^{*}\right)^{*}+\mathbf{T}\mathbf{a}_{n}$

where it was assumed that scattering, as well as transition matrices, are
symmetric and where only formulations involving excitation vectors
$\mathbf{a}_{n}$ are shown for brevity.

Algorithm 1 Matrix-free algorithm

1: $m=0$

2: $\mathbf{a}_{\mathrm{b}}\leftarrow\mathrm{rand}$

3:while stopping criteria are not met do

4: $\mathbf{a}_{m}\leftarrow\mathbf{a}_{m}/|\mathbf{a}_{m}|$

5: modified Gram-Schmidt procedure over $\\{\mathbf{a}_{m}\\}$

6: $\mathbf{f}_{1m}\leftarrow\mathbf{M}\mathbf{a}_{m}$

7:
$\widehat{\mathbf{a}}\leftarrow\mathbf{P}\left(\mathbf{a}_{m}+2\mathbf{f}_{1m}\right)^{*}$

8:
$\widehat{\mathbf{f}}_{2m}\leftarrow\mathbf{M}_{\mathrm{b}}\widehat{\mathbf{a}}$

9:
$\mathbf{f}_{m}\leftarrow\mathbf{f}_{1m}+\mathbf{P}\widehat{\mathbf{f}}_{2m}^{\ast}$

10: $\mathbf{A}_{m}\leftarrow\sum_{p\leq
m}\mathbf{f}_{p}\mathbf{a}_{p}^{\mathrm{H}}$

11: $\\{t_{n}\\}_{m}\leftarrow\mathrm{eig}(\mathbf{A}_{m})$

12: $\mathbf{P}_{m}\leftarrow\sum_{p\leq
m}\mathbf{a}_{p}\mathbf{a}_{p}^{\mathrm{H}}$

13: $\mathbf{a}_{m+1}\leftarrow\mathbf{f}_{m}-\mathbf{P}_{m}\mathbf{f}_{m}$

14: $m\leftarrow m+1$

15:end while

An example of the procedure used to estimate substructure characteristic modes
using transition matrices in matrix-free manner is sketched in Algorithm 1,
where $\mathbf{M}=\mathbf{T}$,
$\mathbf{M}_{\mathrm{b}}=\mathbf{T}_{\mathrm{b}}$,
$\mathbf{A}=2\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}\mathbf{T}+\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T}$
abbreviates the matrix to be decomposed, and $\mathbf{P}=\mathbf{I}$ is an
identity matrix. The algorithm is stopped when magnitude
$\left|\mathbf{a}_{m+1}\right|$ is sufficiently small or relative changes in
estimated eigenvalues $t_{n}$ are sufficiently small. Steps no. 5 and 6 are
the solutions to scattering problems involving full structure and background,
respectively, and can be obtained from any full-wave electromagnetic solver.
The modified Gram–Schmidt procedure is used in Algorithm 1 to assure its
stability [38].

Another possibility to evaluate substructure characteristic modes is to employ
scattering dyadic matrices, accessible with arbitrary electromagnetic solver
[7]. In this case, the matrices $\tilde{\mathbf{S}}$ and
$\tilde{\mathbf{S}}_{\mathrm{b}}$ are defined as in [7, eq. (21)] and denoted
here as scattering and background scattering dyadic matrices, respectively.
The matrices are not transposed symmetric [7, eq. (4)], so the Algorithm 1
must be modified by settings $\mathbf{M}=\tilde{\mathbf{S}}$,
$\mathbf{M}_{\mathrm{b}}=\tilde{\mathbf{S}}_{\mathrm{b}}$, and with matrix
$\mathbf{P}$ being indexing matrix containing zeros except positions $P_{pq}$
and $P_{pq}$ where pairs of the quadrature points are mapped as

$P_{pq}=P_{pq}=1\Longleftrightarrow\hat{\boldsymbol{r}}_{p}=-\hat{\boldsymbol{r}}_{q}.$
(18)

In addition, matrix $\mathbf{P}$ is further modified by flipping the $\pm$
sign for respective positions where the quadrature points in
$\vartheta$-polarization block of the dyadic lie on $\pm z$ axis, and for all
entries corresponding to $\varphi$-polarization block except of points lying
on $\pm z$ axis.

##  Appendix B Equivalence Between MoM-based and Scattering-Based Substructure
Characteristic Modes

The modified transition matrix (14) resembles the transition matrix-based
expression in [5] with the difference that the spherical wave matrix
$\widetilde{\mathbf{U}}$ is complex-valued. Rewriting (14) using block
matrices (10) produces

$\widetilde{\mathbf{T}}=\mathbf{U}_{1}\begin{bmatrix}-\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}\widetilde{\mathbf{Z}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}^{\mathrm{H}}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-\mathrm{H}}&\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}\widetilde{\mathbf{Z}}^{-1}\\\\[5.0pt]
\widetilde{\mathbf{Z}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}^{\mathrm{H}}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-\mathrm{H}}&-\widetilde{\mathbf{Z}}^{-1}\end{bmatrix}\mathbf{U}_{1}^{\mathrm{T}}$
(19)

which partly resembles a block inversion of the MoM matrix (10) except for
some Hermitian transposes.

To express the transition matrices of the composite object $\varOmega$ and
background object $\varOmega_{\mathrm{b}}$ in MoM system matrices, we use [5]

$\mathbf{T}=-\mathbf{U}_{1}\mathbf{Z}^{-1}\mathbf{U}_{1}^{\mathrm{T}}\quad\text{and
}\mathbf{T}_{\mathrm{b}}=-\mathbf{U}_{1}\begin{bmatrix}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}\mathbf{U}_{1}^{\mathrm{T}}.$ (20)

Substituting these T-matrices into (7) and using block matrix inversion
together with algebraic manipulations outlined below we realize that (11),
(13), (7) are all identical.

The derivation starts with reformulation of (7) in MoM matrices

$2\mathbf{T}\mathbf{T}^{\mathrm{H}}_{\mathrm{b}}+\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T}\\\
=\mathbf{U}_{1}\mathbf{Z}^{-1}\left(2\mathbf{R}\begin{bmatrix}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-\mathrm{H}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}-\mathbf{Z}\begin{bmatrix}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-\mathrm{H}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}-\mathbf{1}\right)\mathbf{U}_{1}^{\mathrm{T}},$
(21)

where $\mathbf{R}=\mathbf{U}_{1}^{\mathrm{T}}\mathbf{U}_{1}$ has been used
[26]. Realizing further that $2\mathbf{R}=\mathbf{Z}+\mathbf{Z}^{\mathrm{H}}$
for lossless scatterers, the relation simplifies to

$2\mathbf{T}\mathbf{T}^{\mathrm{H}}_{\mathrm{b}}+\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T}=\mathbf{U}_{1}\mathbf{Z}^{-1}\left(\mathbf{Z}^{\mathrm{H}}\begin{bmatrix}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-\mathrm{H}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{bmatrix}-\mathbf{1}\right)\mathbf{U}_{1}^{\mathrm{T}}\\\
=\mathbf{U}_{1}\mathbf{Z}^{-1}\begin{bmatrix}\mathbf{0}&\mathbf{0}\\\
\mathbf{Z}_{\mathrm{b}\mathrm{c}}^{\mathrm{H}}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-\mathrm{H}}&-\mathbf{1}\end{bmatrix}\mathbf{U}_{1}^{\mathrm{T}}.$
(22)

The final step is the use of block matrix inversion

$\mathbf{Z}^{-1}=\begin{bmatrix}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}+\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}\widetilde{\mathbf{Z}}^{-1}\mathbf{Z}_{\mathrm{c}\mathrm{b}}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}&-\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}\mathbf{Z}_{\mathrm{b}\mathrm{c}}\widetilde{\mathbf{Z}}^{-1}\\\\[5.0pt]
-\widetilde{\mathbf{Z}}^{-1}\mathbf{Z}_{\mathrm{c}\mathrm{b}}\mathbf{Z}_{\mathrm{b}\mathrm{b}}^{-1}&\widetilde{\mathbf{Z}}^{-1}\end{bmatrix}$
(23)

which is identical to $\widetilde{\mathbf{T}}$ in (19)

$2\mathbf{T}\mathbf{T}^{\mathrm{H}}_{\mathrm{b}}+\mathbf{T}_{\mathrm{b}}^{\mathrm{H}}+\mathbf{T}=\widetilde{\mathbf{T}}$
(24)

proving that (11), (13), (7) are the same.

Notice that the above derivation demanded lossless scatterer. For lossy cases,
matrix $\widetilde{\mathbf{R}}$ differs from the radiation operator
$\widetilde{\mathbf{U}}_{1}^{\mathrm{H}}\widetilde{\mathbf{U}}_{1}$. For the
special case with lossless background, the scattering-based formulation of
characteristic modes is identical to the $\widetilde{\mathbf{Z}}$-formulation
on the right in (11).

## References

  * [1] C. G. Montgomery, R. H. Dicke, and E. M. Purcell, _Principles of Microwave Circuits_. New York, NY: McGraw-Hill, 1948\.
  * [2] R. Garbacz, “Modal expansions for resonance scattering phenomena,” _Proc. IEEE_ , vol. 53, no. 8, pp. 856–864, Aug. 1965.
  * [3] R. F. Harrington and J. R. Mautz, “Theory of characteristic modes for conducting bodies,” _IEEE Trans. Antennas Propag._ , vol. 19, no. 5, pp. 622–628, 1971.
  * [4] B. K. Lau, M. Capek, and A. M. Hassan, “Characteristic modes: Progress, overview, and emerging topics,” _IEEE Antennas and Propagation Magazine_ , vol. 64, no. 2, pp. 14–22, 2022.
  * [5] M. Gustafsson, L. Jelinek, K. Schab, and M. Capek, “Unified theory of characteristic modes: Part I–Fundamentals,” _IEEE Trans. Antennas Propag._ , vol. 70, no. 12, pp. 11 801–11 813, 2022.
  * [6] ——, “Unified theory of characteristic modes: Part II–Tracking, losses, and FEM evaluation,” _IEEE Trans. Antennas Propag._ , vol. 70, no. 12, pp. 11 814–11 824, 2022.
  * [7] M. Capek, J. Lundgren, M. Gustafsson, K. Schab, and L. Jelinek, “Characteristic mode decomposition using the scattering dyadic in arbitrary full-wave solvers,” _IEEE Trans. Antennas Propag._ , vol. 71, no. 1, pp. 830–839, 2023.
  * [8] J. Lundgren, K. Schab, M. Capek, M. Gustafsson, and L. Jelinek, “Iterative calculation of characteristic modes using arbitrary full-wave solvers,” _IEEE Antennas Wireless Propag. Lett._ , vol. 22, no. 4, pp. 799–803, 2023\.
  * [9] J. Ethier and D. McNamara, “Sub-structure characteristic mode concept for antenna shape synthesis,” _Electronics letters_ , vol. 48, no. 9, p. 1, 2012\.
  * [10] H. Alroughani, “An appraisal of the characteristic modes of composite objects,” Ph.D. dissertation, University of Ottawa (Canada), 2013.
  * [11] G. Angiulli and G. Di Massa, “Scattering from arbitrarily shaped microstrip patch antennas using the theory of characteristic modes,” in _IEEE Antennas and Propagation Society International Symposium and USNC/URSI National Radio Science Meeting_ , ser. APS-98. IEEE.
  * [12] R. Zhao, Y. Lu, G. S. Cheng, W. Zhu, J. Hu, and H. Bagci, “Sub-structure characteristic mode analysis of microstrip antennas using a global multitrace formulation,” _IEEE Trans. Antennas Propag._ , vol. 71, no. 12, pp. 10 026–10 031, 2023.
  * [13] R. Luomaniemi, P. Ylä-Oijala, A. Lehtovuori, and V. Viikari, “Designing hand-immune handset antennas with adaptive excitation and characteristic modes,” _IEEE Trans. Antennas Propag._ , vol. 69, no. 7, pp. 3829–3839, 2021\.
  * [14] P. Parhami, Y. Rahmat-Samii, and R. Mittra, “Technique for calculating the radiation and scattering characteristics of antennas mounted on a finite ground plane,” in _Proceedings of the Institution of Electrical Engineers_ , vol. 124, no. 11. IET, 1977, pp. 1009–1016.
  * [15] T. Cwik, J. Patterson, and T. Lockhart, “Constructing matrix Green’s functions for radiation and scattering problems,” in _Digest on Antennas and Propagation Society International Symposium_ , 1989, pp. 586–589 vol.2.
  * [16] Q. I. Dai, H. Gan, C. C. Weng, and C.-F. Wang, “Characteristic mode analysis using Green’s function of arbitrary background,” in _2016 IEEE International Symposium on Antennas and Propagation (APSURSI)_. IEEE, 2016, pp. 423–424.
  * [17] A. Alakhras and D. McNamara, “Sub-structure characteristic mode computation utilising field-based MM/GTD hybrid methods,” _Journal of Electromagnetic Waves and Applications_ , vol. 34, no. 13, pp. 1812–1821, 2020\.
  * [18] R. F. Harrington, _Field Computation by Moment Methods_. New York, NY: Macmillan, 1968.
  * [19] H. Alroughani, J. L. T. Ethier, and D. A. McNamara, “On the classification of characteristic modes, and the extension of sub-structure modes to include penetrable material,” in _2014 International Conference on Electromagnetics in Advanced Applications (ICEAA)_. IEEE, 2014.
  * [20] G. Kristensson, _Scattering of Electromagnetic Waves by Obstacles_. Edison, NJ: SciTech Publishing, an imprint of the IET, 2016.
  * [21] D. M. Pozar, _Microwave Engineering_ , 3rd ed. New York, NY: John Wiley & Sons, 2005.
  * [22] M. I. Mishchenko, L. D. Travis, and D. W. Mackowski, “T-matrix computations of light scattering by nonspherical particles: A review,” _J. Quant. Spectrosc. Radiat. Transfer_ , vol. 55, no. 5, pp. 535–575, 1996.
  * [23] M. Capek and K. Schab, “Computational aspects of characteristic mode decomposition: An overview,” _IEEE Antennas and Propagation Magazine_ , vol. 64, no. 2, pp. 23–31, 2022.
  * [24] R. J. Garbacz and R. H. Turpin, “A generalized expansion for radiated and scattered fields,” _IEEE Trans. Antennas Propag._ , vol. 19, no. 3, pp. 348–358, 1971.
  * [25] K. Schab, F. W. Chen, L. Jelinek, M. Capek, J. Lundgren, and M. Gustafsson, “Characteristic modes of frequency-selective surfaces and metasurfaces from S-parameter data,” _IEEE Trans. Antennas Propag._ , vol. 71, no. 12, pp. 9696–9706, 2023.
  * [26] D. Tayli, M. Capek, L. Akrou, V. Losenicky, L. Jelinek, and M. Gustafsson, “Accurate and efficient evaluation of characteristic modes,” _IEEE Trans. Antennas Propag._ , vol. 66, no. 12, pp. 7066–7075, 2018.
  * [27] V. Losenicky, L. Jelinek, M. Capek, and M. Gustafsson, “Method of moments and T-matrix hybrid,” _IEEE Trans. Antennas Propag._ , vol. 70, no. 5, pp. 3560–3574, 2022.
  * [28] L. Jelinek, M. Gustafsson, K. Schab, M. Capek, and E. Moreno, “Transition matrix in characteristic modes theory,” in _2022 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting (AP-S/URSI)_. IEEE, jul 2022.
  * [29] M. Masek, M. Capek, L. Jelinek, and K. Schab, “Modal tracking based on group theory,” _IEEE Trans. Antennas Propag._ , vol. 68, no. 2, pp. 927–937, Feb. 2020.
  * [30] (2022) Altair FEKO. Altair. [Online]. Available: https://www.altair.com/feko
  * [31] M. Capek. (2024) Scattering dyadic characteristic modes (with Altair FEKO). [Online]. Available: https://github.com/kschab/scattering-dyadic-characteristic-modes/tree/main/FEM_MoM-FEKO
  * [32] (2017) Cubesat 101 basic concepts and processes for first-time cubesat developers. National Aeronautics and Space Administration (NASA). [Online]. Available: https://www.nasa.gov/wp-content/uploads/2017/03/nasa_csli_cubesat_101_508.pdf
  * [33] X. Deng, Y. Chen, and S. Yang, “Characteristic mode formulation for antennas with waveguide port feeding structures,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 20, no. 10, pp. 2063–2067, 2021.
  * [34] W. Kahn and H. Kurss, “Minimum-scattering antennas,” _IEEE Trans. Antennas Propag._ , vol. 13, no. 5, pp. 671–675, 1965.
  * [35] P. Rogers, “Application of the minimum scattering antenna theory to mismatched antennas,” _IEEE Trans. Antennas Propag._ , vol. 34, no. 10, pp. 1223–1228, 1986.
  * [36] Y. G. Kim and S. Nam, “Determination of the generalized scattering matrix of an antenna from characteristic modes,” _IEEE Trans. Antennas Propag._ , vol. 61, no. 9, pp. 4848–4852, 2013.
  * [37] M. Alian and N. Noori, “A domain decomposition method for the analysis of mutual interactions between antenna and arbitrary scatterer using generalized scattering matrix and translation addition theorem of SWFs,” _IEEE Trans. Antennas Propag._ , vol. 71, no. 10, p. 8088–8096, Oct. 2023.
  * [38] G. H. Golub and C. F. van Loan, _Matrix Computations_ , 4th ed. Baltimore, MD: The Johns Hopkins University Press, 2013.

|  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.
---|---
|  Lukas Jelinek was born in Czech Republic in 1980. He 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.
---|---
|  Miloslav Capek (M’14, SM’17) received the M.Sc. degree in Electrical
Engineering 2009, the Ph.D. degree in 2014, and was appointed a Full Professor
in 2023, 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 include electromagnetic theory, electrically small
antennas, antenna design, numerical techniques, and optimization. He authored
or co-authored over 160 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. He received the IEEE Antennas and
Propagation Edward E. Altshuler Prize Paper Award 2023.
---|---
|  Johan Lundgren (M’22) is a postdoctoral researcher at Lund University. he
received his M.Sc degree in engineering physics 2016 and Ph.D. degree in
Electromagnetic Theory in 2021 all from Lund University, Sweden. His research
interests are in electromagnetic scattering, wave propagation, computational
electromagnetics, functional structures, meta-materials, inverse scattering
problems, imaging, and measurement techniques.
---|---
|  Kurt Schab (S’09, M’16) is an Assistant Professor of Electrical Engineering
at Santa Clara University, Santa Clara, CA USA. He received the B.S. degree in
electrical engineering and physics from Portland State University in 2011 and
the M.S. and Ph.D. degrees in electrical engineering from the University of
Illinois at Urbana-Champaign in 2013 and 2016, respectively. From 2016 to 2018
he was a Postdoctoral Research Scholar at North Carolina State University in
Raleigh, North Carolina. His research focuses on the intersection of numerical
methods, electromagnetic theory, and antenna design.
---|---