# Circle Packings from Tilings of the Plane

Philip Rehwinkel, Ian Whitehead, David Yang, Mengyuan Yang

###### Abstract.

We introduce a new class of fractal circle packings in the plane, generalizing
the polyhedral packings defined by Kontorovich and Nakamura. The existence and
uniqueness of these packings are guaranteed by infinite versions of the Koebe-
Andreev-Thurston theorem. We prove structure theorems giving a complete
description of the symmetry groups for these packings. And we give several
examples to illustrate their number-theoretic and group-theoretic
significance.

##  1\. Introduction

The well-known Apollonian circle packing can be constructed from a set of four
base circles, and four dual circles, as shown in Figure 1. The orbit of the
base circles under the group generated by reflections through the dual circles
is the packing, an infinite fractal set of circles. Beyond their aesthetic
appeal, Apollonian packings have properties of great interest in number
theory, group theory, and fractal geometry.

(a) Base (blue) and dual (red) circles

(b) Apollonian packing

Figure 1. Constructing the Apollonian packing

In [13], Kontorovich and Nakamura define polyhedral circle packings,
generalizing the Apollonian packing construction. Any circle configuration has
a tangency graph, with a vertex for each circle and an edge for each tangency
between circles. In the case of the Apollonian packing, both the base circles
and the dual circles have tetrahedral tangency graphs. In general, one can
start with a finite set of base circles whose tangency graph is the graph of
any polyhedron and a finite set of dual circles whose tangency graph is the
graph of the dual polyhedron. The orbit of the base circles under the group
generated by the dual circles is a polyhedral packing. Polyhedral packings
encompass many of the generalizations of the Apollonian packing that have been
studied previously. For example, the packing introduced by Guettler and
Mallows [9] is the octahedral packing; the $\mathbb{Q}[\sqrt{-2}]$ packing
studied by Stange [22] is the cubic packing.

In this article, we study packings which originate from infinite
configurations of base and dual circles, a further generalization. A
particularly symmetric example is shown in Figure 2. In this example, both the
base and dual circles have the square lattice as their tangency graphs. Again,
the orbit of the base circles under the group generated by the dual circles is
a fractal set of circles. We call this object the square packing.

(a) Base (blue) and dual (red) circles

(b) Square packing

Figure 2. Constructing the square packing

In general, we work with a base circle configuration $B$ and a dual
configuration $\hat{B}$ whose tangency graphs determine dual tilings, or
cellular decompositions of the sphere $\hat{\mathbb{C}}$ or the plane
$\mathbb{C}$–see Definition 2.1. Tilings of the sphere give rise to polyhedral
packings, while tilings of the plane give rise to new examples. In both cases,
the final packings have similar geometric properties. For example, the circles
in the packing are pairwise disjoint or tangent, they can be oriented with
disjoint interiors, and the interiors are dense in the ambient space. The
symmetry groups of polyhedral packings and of our new examples have similar
structure, as we prove in Section 4. And some of our new examples have
integrality properties, raising number-theoretic questions.

One motivation to introduce these new packings comes from the literature on
the Koebe-Andreev-Thurston theorem and its generalizations. The Koebe-Andreev-
Thurston theorem [12, 3, 24] is the remarkable result that any pair of finite
graphs $G$, $\hat{G}$ representing a polyhedron and its dual can be realized
as the tangency graphs for a pair of dual circle configurations $B$,
$\hat{B}$. Moreover, the circle configurations $B$, $\hat{B}$ are unique up to
conformal automorphism of $\hat{\mathbb{C}}$ (see Theorem 3.1 for a precise
statement). The Koebe-Andreev-Thurston theorem implies the existence and
uniqueness up to Möbius transformation of a circle packing for every
polyhedron.

It is natural to try to extend this theorem to infinite graphs. Important work
of Beardon-Stephenson and of Schramm achieves this in many cases [4, 19, 20],
following a constructive approach suggested by Thurston. Stephenson’s text
[23] synthesizes this work. One result is that any infinite graph $G$
representing a triangulation of the plane can be realized as the tangency
graph for a circle configuration $B$, and $B$ is unique up to conformal
automorphism of $\mathbb{C}$ (see Theorem 3.2). The proofs use deep geometric
ideas—mappings between circle configurations give a discrete analogue of the
Riemann mapping theorem.

Just as the construction of polyhedral circle packings relies on the finite
Koebe-Andreev-Thurston theorem, our construction relies on its infinite
generalizations. In Section 3, we state various versions of the theorem which
imply the existence and uniqueness up to conformal automorphism of many of our
packings. We also make a more general conjecture which would imply existence
and uniqueness in all cases. Our work combines geometric ideas from the
infinite Koebe-Andreev-Thurston theorem with arithmetic ideas from the
Apollonian packing and its relatives.

Another motivation comes from the definitions of crystallographic packings in
[13]. These are a class of packings which generalize the Apollonian packing,
and encompass many known examples of circle and sphere packings. To construct
crystallographic circle packings, one can start with a geometrically finite
reflection group which acts on $\mathbb{H}^{3}$ with finite covolume. Each
wall of the fundamental chamber intersects the spherical boundary of
$\mathbb{H}^{3}$ in a circle. Suppose that these circles are partitioned into
two sets, a “cluster” and “cocluster,” such that circles in the cluster are
pairwise disjoint or tangent, and each circle in the cluster is disjoint,
tangent, or orthogonal to each circle in the cocluster. Then the orbit of the
cluster under the group generated by reflections across the cocluster is a
crystallographic packing. All polyhedral packings are crystallographic, but
not all crystallographic packings are polyhedral.

Kontorovich and Nakamura classify superintegral crystallographic circle
packings. Every superintegral packing arises from an arithmetic finite-
covolume reflection group acting on $\mathbb{H}^{3}$. There are finitely many
such groups up to commensurability, tabulated in [18], and all the non-
cocompact groups in the tabulation give rise to packings. The most interesting
case is the Bianchi group of the Eisenstein integers. This is represented by
the Coxeter-Dynkin diagram:

Each vertex represents a circle, and the edge types indicate angles between
circles. There is a finite-index subgroup with the diagram:

See Figure 9 (A), where this diagram is realized as a set of five circles. If
the left vertex is the cluster and the remaining vertices are the cocluster,
the result is the triangular packing of Figure 10. Similarly, if the right
vertex is the cluster and the remaining vertices are the cocluster, the result
is the hexagonal packing of Figure 11. These are fundamental examples of
superintegral crystallographic packings. We conjecture that they are not
polyhedral, or commensurate (on the level of hyperbolic reflection groups) to
any polyhedral packing. Nevertheless, our construction allows us to work with
these packings similarly to polyhedral packings. One might hope to realize all
superintegral crystallographic circle packings with our construction.

A final motivation comes from limits in polyhedral packing families. In [2],
fractal dimensions for many polyhedral packings are computed. It is observed
that for some sequences of polyhedra–pyramids, prisms, antiprisms, etc.–the
fractal dimensions converge to a limit. More surprisingly, the packings
themselves converge to a well-defined limit. This phenomenon is illustrated in
Figure 3. A configuration of circles in the polyhedral packing for the
100-sided prism is shown. As the number of sides increases, this configuration
approaches a configuration found in the square packing of Figure 2. The limit
of the prism packings is our square packing.

(a) 100-sided prism configuration

(b) Square configuration

Figure 3. Limiting configurations

The article [2] focuses on examples, and does not give a general criterion for
a sequence of circle packings to converge. Some of the limit packings are also
polyhedral–e.g. the limit of the pyramid packings is the original Apollonian
packing–and others are not polyhedral but satisfy our definition. We discuss
all these examples in Section 5. It would be very interesting to know whether
our definition gives the closure of the space of polyhedral packings. This
would require a more systematic development of the notion of a limit in this
space.

This article is structured as follows. In Section 2, Definitions 2.1 and 2.7
describe the class of packings we study. We prove some of the fundamental
geometric consequences of these definitions, and compare them to the
definitions of polyhedral, crystallographic, and Kleinian packings. In Section
3, we recall versions of the Koebe-Andreev-Thurston theorem and propose our
extension, Conjecture 3.3. In Section 4, we analyze the symmetries of our
packings. Theorems 4.1, 4.2, and 4.10 give a complete description of the
symmetry group. In Section 5, we give three main examples of our packings,
which we call the triangular, square, and hexagonal packings. We focus on
their arithmetic properties–integrality, quadratic and linear forms. In
Section 6, we give a broader class of examples, with a focus on symmetries.
Theorem 6.2 shows that all 17 wallpaper groups appear in the symmetry groups
of packings.

In future work, we hope to elaborate on the number theory of the packings
described here. As indicated in Section 5, many of our packings have integral
curvatures. Is there an asymptotic formula for curvatures in the packing, as
in [14]? Is there a local-to-global principle for curvatures, as conjectured
in [8]? These questions can be answered by the methods of [7] for some
periodic packings. They become much more subtle if the base configuration
lacks any symmetry. It is possible to construct integral, aperiodic packings
using the same refinement method as in the proof of Theorem 6.2. These
packings are geometrically interesting, and their number theory remains to be
explored.

###  1.1. Acknowledgements

We thank Daniel Allcock, Arthur Baragar, Michael Dougherty, Cathy Hsu, Anna
Felickson, Alex Kontorovich, Alice Mark, and Kate Stange for helpful
conversations relating to this project. We thank Nooria Ahmed, William Ball,
Ellis Buckminster, Emilie Rivkin, Dylan Torrance, Jake Viscusi, Runze Wang,
and Gary Yang for raising some of the questions that influenced our work. We
are grateful to Swarthmore College for funding the summer research project
that led to this article, and we thank all the faculty and students in the
Department of Mathematics and Statistics who helped create a productive
research community.

##  2\. Definitions and Basic Properties

The extended complex plane $\hat{\mathbb{C}}$ is $\mathbb{C}\cup\\{\infty\\}$,
with the topology of the sphere. Our packings consist of oriented generalized
circles, i.e. circles and lines, in $\hat{\mathbb{C}}$. Each generalized
circle divides $\hat{\mathbb{C}}$ into two simply connected regions. Choosing
an orientation for the circle is equivalent to choosing one of these regions
to be the interior, and the other to be the exterior. When we refer to circles
in this article, we always mean oriented generalized circles. Oriented
generalized circles have an action by the group of holomorphic and
antiholomorphic Möbius transfomations
$\mathrm{M\ddot{o}b}\cong\mathrm{SL}_{2}(\mathbb{C})\rtimes\mathbb{Z}/2\mathbb{Z}$.
For more details on this setup, see [15].

For any collection of circles $B$, we can associate its tangency graph $G_{B}$
which has a vertex for every circle and an edge between each pair of tangent
circles. We will say that a collection of circles accumulates at a point
$x\in\hat{\mathbb{C}}$ if any open neighborhood of $x$ contains infinitely
many circles from this collection. We now define the notions of base and dual
circle configurations, which are the starting point for the packings we
construct.

######  Definition 2.1.

Let $B$ and $\hat{B}$ be two collections of oriented generalized circles, with
tangency graphs $G_{B}$ and $G_{\hat{B}}$, respectively. Then $B$ is called a
base configuration and $\hat{B}$ is called a dual configuration if the
following properties hold:

  1. (1)

The circles in $B$ are pairwise disjoint or tangent, with disjoint interiors,
and the same holds for the circles in $\hat{B}$.

  2. (2)

The tangency graphs $G_{B}$ and $G_{\hat{B}}$ are each nontrivial, connected,
and are duals of each other.

  3. (3)

If a circle in $B$ and a circle $\hat{B}$ intersect, they do so orthogonally
and they correspond to a face-vertex pair in the tangency graphs. Otherwise,
their interiors are disjoint.

  4. (4)

$B\cup\hat{B}$ has at most one accumulation point.

Note that the roles of $B$ and $\hat{B}$ are interchangeable in this
definition. An example of a base and dual configuration pair is shown in
Figure 4.

Figure 4. Base and dual circle configurations

The tangency graphs $G_{B}$ and $G_{\hat{B}}$ are simple, and by (1) they have
natural planar embeddings, placing each vertex at the center of its circle (or
choosing an arbitrary interior point as the vertex if the circle is a line).
The duality between graphs, with faces of $G_{B}$ corresponding to vertices of
$G_{\hat{B}}$ and vice versa, is defined in this context. The four properties,
and especially this duality, imply stronger geometric statements, for example:

######  Proposition 2.2.

Suppose that two circles $c_{1},c_{2}\in B$ are tangent at a point
$x\in\hat{\mathbb{C}}$. Then no other circles from $B$ go through this point.
Moreover, if $x$ is not the accumulation point, then it is also a point of
tangency for exactly two dual circles $d_{1},d_{2}\in\hat{B}$, which intersect
$c_{1},c_{2}$ orthogonally.

###### Proof.

First, there cannot be any additional circle in $B$ tangent to $c_{1},c_{2}$
at $x$ as this would violate the disjoint interiors property (1). Next,
suppose that $x$ is not the accumulation point, and consider the edge
connecting $c_{1},c_{2}$ in $G_{B}$. This corresponds to an edge connecting
some pair of tangent dual circles $d_{1},d_{2}\in G_{\hat{B}}$, which
intersect $c_{1},c_{2}$ orthogonally. After applying a Möbius transformation,
we may assume that $c_{1}$ and $c_{2}$ are parallel horizontal lines which are
tangent at the point $x=\infty$. Then it is clear that $d_{1},d_{2}$ must be
parallel vertical lines which are tangent at the same point $x$. Again, there
cannot be any additional circle in $\hat{B}$ tangent to $d_{1},d_{2}$ at $x$
as this would violate property (1). ∎

######  Proposition 2.3.

Each circle in $c\in B$ is orthogonal to at least three dual circles in
$\hat{B}$. These dual circles can be labeled by elements of $\mathbb{Z}$ or
$\mathbb{Z}/n\mathbb{Z}$ so that consecutive circles are tangent.

###### Proof.

The vertex corresponding to $c$ in $G_{B}$ must be connected to the rest of
the graph by at least three edges; if it were connected by two edges, the dual
graph would have a double edge, and if it were connected by one edge, then the
dual graph would have a loop. Thus the vertex is incident to at least three
faces. These faces must correspond to distinct dual circles because each dual
circle is uniquely determined by two points of intersection with $c$. Then the
rest of the statement follows from the fact that the faces incident to a
vertex in a planar graph can be ordered cyclically so that any two consecutive
faces are adjacent. ∎

We will describe the situation of Proposition 2.3 by saying that $c$ is ringed
by circles from $\hat{B}$. One subtlety in the proof is the possibility that
$c$ could be ringed by infinitely many circles. This can happen if the unique
accumulation point lies on $c$, and in this case, the orthogonal circles can
be labeled by the integers so that any two consecutive circles are tangent.

######  Proposition 2.4.

The circles in $B$, $\hat{B}$, and their interiors cover all points in
$\hat{\mathbb{C}}$ other than the accumulation point.

The accumulation point may or may not be covered, depending on whether or not
it lies on a circle.

###### Proof.

It suffices to show that all points on one face in the planar embedding of
$G_{B}$ are covered. After a Möbius transformation, we may assume that this
face does not contain the point $\infty$, and that none of the generalized
circles in $B$ surrounding it are lines. By Proposition 2.3, the face is a
simply connected polygon (possibly with infinitely many sides), with vertices
at the centers of circles in $B$. Each edge goes through a point of tangency
between two circles in $B$, and is orthogonal to both circles. Then by
Proposition 2.2, the dual circle in $\hat{B}$ corresponding to the face is
tangent to each edge, i.e. it is inscribed in the polygon. Removing the dual
circle and its interior leaves one connected component for each vertex, and
each connected component is contained in the interior of the corresponding
circle. ∎

A graph is is said to be $n$-connected if after removing any $n-1$ vertices
and their adjacent edges, the graph remains connected. It is said to be
$n$-edge connected if after removing any $n-1$ edges, the graph remains
connected. Because $G_{B}$ and $G_{\hat{B}}$ are a pair of dual simple planar
graphs, they are necessarily 3-edge connected. If removing an edge could
disconnect $G_{B}$, then $G_{\hat{B}}$ would have to contain a loop, and if
removing two edges could disconnect $G_{B}$, then $G_{\hat{B}}$ would have to
contain a double edge. In fact, we can make an even stronger connectedness
statement:

######  Proposition 2.5.

The tangency graphs $G_{B}$ and $G_{\hat{B}}$ are $3$-connected.

###### Proof.

It suffices to show that if one or two vertices and their adjacent edges are
removed from $G_{B}$, then all faces in the planar embedding of $G_{B}$ remain
simply connected. This implies that $G_{B}$ remains connected.

If a single vertex $c_{1}$ is removed, then the faces incident to that vertex
are all identified. By Proposition 2.3, these faces correspond to the ring of
dual circles around $c_{1}$. Because these faces are all distinct and simply
connected, gluing them at the vertex and along their common edges results in a
larger simply connected face.

If two non-adjacent vertices $c_{1},c_{2}$ are removed, then at each vertex
the incident faces are identified. If $c_{1},c_{2}$ have no face in common,
then the argument is exactly the same as above. If they have a face in common,
there can only be one such face. Indeed, after applying a Möbius
transformation, we may assume that $c_{1},c_{2}$ are concentric circles
centered at the origin. Any dual circle orthogonal to both must be a line
through the origin. By the disjoint interiors property (1), $\hat{B}$ can
contain at most one such line. In this case, removing $c_{1}$ and $c_{2}$
creates one new larger face, the union of all faces incident to $c_{1}$ or
$c_{2}$. Because the union of faces incident to $c_{1}$ is simply connected,
the union of faces incident to $c_{2}$ is simply connected, and they overlap
in a unique simply connected face, the new face must be simply connected.

Finally, if two adjacent vertices $c_{1},c_{2}$ are removed, by Proposition
2.2, there are exactly two faces incident to these two vertices, meeting along
the edge from $c_{1}$ to $c_{2}$. The union of these two faces is simply
connected. Removing $c_{1}$ and $c_{2}$ creates one new larger face, the union
of all faces incident to $c_{1}$ or $c_{2}$. Because the union of faces
incident to $c_{1}$ is simply connected, the union of faces incident to
$c_{2}$ is simply connected, and they overlap in a simply connected union of
two faces, the new face must be simply connected. ∎

By the compactness of $\hat{\mathbb{C}}$, $B$ is finite if an only if it has
no accumulation point. In this case $G_{B}$ is a finite 3-connected simple
planar graph, so by Steinitz’s theorem, it is the graph of a polyhedron.

One further property we might ask for in base and dual configurations is
periodicity:

######  Definition 2.6.

For $n=1,\,2$, $B$ is $n$-periodic, i.e. periodic under an $n$-dimensional
lattice, if there exist $v_{1},\dots,v_{n}\in\mathbb{C}$, linearly independent
over $\mathbb{R}$, such that $B+v_{1}=\dots=B+v_{n}=B$.

If $B$ is $1$-periodic or $2$-periodic, it necessarily has an accumulation
point at $\infty$. The tangency graph of a $2$-periodic $B$ is a $2$-periodic
tiling of the plane.

We are now ready to define the class of circle packings that we will study.
For any circle $d$, let $\sigma_{d}$ denote the reflection across $d$, a
Möbius transformation.

######  Definition 2.7.

The packing $\mathscr{P}$ is the orbit of $B$ under the group generated by
reflections $\sigma_{d}$ across circles $d\in\hat{B}$. The dual packing
$\hat{\mathscr{P}}$ is the orbit of $\hat{B}$ under the same group.

The superpacking is the orbit of $B$ under the group generated by reflections
across circles in $B$ and $\hat{B}.$

Note that under this definition, the dual packing and superpacking contain two
oppositely-oriented copies of each circle. It is possible to assign a single
orientation to circles in the dual packing and the superpacking in a
consistent way, but since we will not use the orientations of these circles in
an important way, we do not pursue this. We will only be concerned with the
orientations of the circles in $\mathscr{P}$.

For each packing $\mathscr{P}$, we define the following symmetry groups:

  1. (1)

$\Gamma=\textsc{Sym}(\mathscr{P},\hat{\mathscr{P}})$: the group of Möbius
transformations that preserve both the packing and the dual packing;

  2. (2)

$\Gamma_{1}=\langle\sigma_{d}\,:\,d\in\hat{B}\rangle$: the group generated by
reflections across the dual circles;

  3. (3)

$\Gamma_{2}=\textsc{Sym}(B,\hat{B})$: the group of Möbius transformations that
preserve both the base configuration and the dual configuration.

Some fundamental geometric properties of $\mathscr{P}$ can be deduced directly
from the definition.

######  Proposition 2.8.

Every circle in $\mathscr{P}$ is either a base circle in $B$ or inside some
dual circle, and every circle in $\hat{\mathscr{P}}$ is either a dual circle
in $\hat{B}$ or inside some dual circle.

Note that when we say circle $c_{1}$ is inside circle $c_{2}$, we only mean
that $c_{1}$ is contained in the union of $c_{2}$ and its interior, not
necessarily that the two interiors are nested.

###### Proof.

By definition, each circle in $\mathscr{P}$ is in $B$ or it is
$\sigma_{d_{1}}\cdots\sigma_{d_{k}}(c)$ for some $c\in B$, $d_{1},\ldots
d_{k}\in\hat{B}$, with consecutive $d_{i}$ distinct. We may assume that $c$ is
not orthogonal to $d_{k}$, because then $d_{k}$ could be dropped from this
expression. Then by property (3), $c$ is outside $d_{k}$, so
$\sigma_{d_{k}}(c)$ is inside $d_{k}$. By property (1), this implies that
$\sigma_{d_{k}}(c)$ is outside $d_{k-1}$, so
$\sigma_{d_{k-1}}\sigma_{d_{k}}(c)$ is inside $d_{k-1}$. Repeating the
argument inductively, we see that $\sigma_{d_{1}}\cdots\sigma_{d_{k}}(d)$ is
inside $d_{1}$. The proof of the second statement is similar. ∎

######  Proposition 2.9.

The circles in $\mathscr{P}$ are pairwise disjoint or tangent, with disjoint
interiors.

###### Proof.

Suppose that we have two circles in $\mathscr{P}$ with overlapping interiors.
After applying a Möbius transformation, we may assume that one of the circles,
$c_{1}$ is in $B$. We may write the other circle as
$\sigma_{d_{1}}\cdots\sigma_{d_{k}}(c_{2})$ for $c_{2}\in B$, $d_{1},\ldots
d_{k}\in\hat{B}$, with consecutive $d_{i}$ distinct. If $d_{k}$ is orthogonal
to $c_{2}$, then we can shorten this expression to
$\sigma_{d_{1}}\cdots\sigma_{d_{k-1}}(c_{2})$. Otherwise, by property (3) of
Definition 2.1, $d_{k}$ and $c_{2}$ have disjoint interiors, and after
inversion, the interior of $\sigma_{d_{k}}(c_{2})$ is contained inside the
interior of $d_{k}$. After the remaining inversions, by property (1), the
interior of $\sigma_{d_{1}}\cdots\sigma_{d_{k}}(c_{2})$ is contained inside
the interior of $d_{1}$. By property (3), in order for $c_{1}$ to intersect
$\sigma_{d_{1}}\cdots\sigma_{d_{k}}(c_{2})$, $c_{1}$ must be orthogonal to
$d_{1}$. Then inverting across $d_{1}$, we find that $c_{1}$ also intersects
$\sigma_{d_{2}}\cdots\sigma_{d_{k}}(c_{2})$. In either case, we have shortened
the string $\sigma_{d_{1}}\cdots\sigma_{d_{k}}$. Repeating this process, we
eventually find two circles $c_{1},c_{2}\in B$ with overlapping interiors,
contradicting property (1). ∎

######  Proposition 2.10.

The interiors of the circles in $\mathscr{P}$ are dense in $\hat{\mathbb{C}}$.

###### Proof.

For $x\in\hat{\mathbb{C}}$, by Proposition 2.4, either $x$ is the accumulation
point, $x$ is in some circle in $B$, or $x$ is in some circle in $\hat{B}$. It
follows that either $x$ is in the closure of the interiors of the circles in
$B$ or $x$ is in the interior of some circle in $\hat{B}$. If $x$ is in the
interior of a dual circle $d_{1}$, reflect it across $d_{1}$. Then either
$\sigma_{d_{1}}(x)$ is in the closure of the interiors of the circles in $B$,
which means that $x$ is in the closure of the interiors of the circles in
$\sigma_{d_{1}}(B)$, or $\sigma_{d_{1}}(x)$ is in the interior of some circle
$d_{2}\in\hat{B}$, which means $x$ is in the interior of
$\sigma_{d_{1}}(d_{2})$. In the latter case, reflect $\sigma_{d_{1}}(x)$
across $d_{2}$ and repeat the process. If, at some step, we find that $x$ is
in the closure of the interiors of the circles in
$\sigma_{d_{1}}\cdots\sigma_{d_{k}}(B)$, then since these are circles in
$\mathscr{P}$, $x$ is in the closure of the interiors of the circles in
$\mathscr{P}$ as desired.

Otherwise, this construction produces an infinite sequence

$\tilde{d}_{1}=d_{1},\,\tilde{d}_{2}=\sigma_{d_{1}}(d_{2}),\,\tilde{d}_{3}=\sigma_{d_{1}}\sigma_{d_{2}}(d_{3}),\ldots$

with consecutive $d_{i}$ distinct, of nested circles in $\hat{\mathscr{P}}$
whose interiors contain $x$. We will show that these circles converge to $x$
in $\hat{\mathbb{C}}$, in the sense that they eventually lie within any open
neighborhood of $x$. After applying a Möbius transformation, we may assume
that $\tilde{d}_{1}$ is not a line and is oriented inward so that its interior
does not contain $\infty$. Since the circles $\tilde{d}_{k}$ are nested, they
all have these properties.

We have that

$\tilde{d}_{k+1}=\sigma_{\tilde{d}_{k}}(\sigma_{d_{1}}\cdots\sigma_{d_{k-1}}(d_{k+1}))$

by Lemma 4.5 below. If $d_{k+1}\neq d_{k-1}$, then
$\sigma_{d_{1}}\cdots\sigma_{d_{k-1}}(d_{k+1})$ and $\tilde{d}_{k}$ are both
contained inside $\tilde{d}_{k-1}$, so neither one contains infinity. Thus,
the reflection of $\sigma_{d_{1}}\cdots\sigma_{d_{k-1}}(d_{k+1})$ across
$\tilde{d}_{k}$ does not contain the center of $\tilde{d}_{k}$, so the radius
of $\tilde{d}_{k+1}$ is at most half the radius of $\tilde{d}_{k}$. Therefore,
if $d_{k+1}\neq d_{k-1}$ for infinitely many values of $k$, then the radii
must approach $0$, so the circles must approach $x$.

If $d_{k+1}=d_{k-1}$, then
$\tilde{d}_{k+1}=\sigma_{\tilde{d}_{k}}(\tilde{d}_{k-1})$ (with the
orientation reversed). If this holds for all but finitely many values of $k$,
then there exists some $K\in\mathbb{N}$ such that for all $k\geq K$,
$\tilde{d}_{k+1}=\sigma_{\tilde{d}_{k}}(\tilde{d}_{k-1})$. Suppose that
$\tilde{d}_{K}$ and $\tilde{d}_{K+1}$ are disjoint. Then, after another Möbius
transformation, we may assume that $\tilde{d}_{K}$ has radius $1$,
$\tilde{d}_{K+1}$ has radius $r<1$, and they are concentric. An inductive
argument shows that the radius of $\tilde{d}_{K+k}$ is $r^{k}$. These radii
approach $0$, so the circles must approach $x$.

On the other hand, suppose that $\tilde{d}_{K}$ and $\tilde{d}_{K+1}$ are
tangent. After a Möbius transformation, we may assume that $\tilde{d}_{K}$ has
radius $1$, and $\tilde{d}_{K+1}$ has radius $r<1$. In this case, an inductive
argument shows that the radius of $\tilde{d}_{K+k}$ is $\frac{1}{1+k(1/r-1)}$.
These radii approach $0$, so the circles must approach $x$.

Since the circles $\tilde{d}_{k}\in\hat{\mathscr{P}}$ approach $x$, there is a
corresponding sequence of circles $\tilde{c}_{k}\in\mathscr{P}$, with
$\tilde{c}_{k}$ orthogonal to $\tilde{d}_{k}$, such that the interiors of the
circles $\tilde{c}_{k}$ come arbitrarily close to $x$. ∎

The residual set of $\mathscr{P}$ is the set of points not in the interior of
any circle.

###  2.1. Relation to other Packing Definitions

In [13] and [10], definitions are given for the related notions of Kleinian,
crystallographic, and polyhedral packings. To compare our definition to these,
we give sufficient conditions for our packings to be polyhedral,
crystallographic, or Kleinian.

The definition of a polyhedral packing coincides with our definition, with the
stricter additional assumption that $G_{B}$ and $G_{\hat{B}}$ are the graphs
of a convex polyhedron and its dual. By the Koebe-Andreev-Thurston theorem,
every polyhedron gives rise to a polyhedral packing. Proposition 2.5 and
Steinitz’s theorem imply the following:

######  Proposition 2.11.

If $B$ and $\hat{B}$ are finite circle configurations satisfying the
conditions of Definition 2.1, then they give rise to a polyhedral circle
packing $\mathscr{P}$.

The definitions of crystallographic and Kleinian packings are more general,
and relate the packing to a discrete group action on a higher-dimensional
hyperbolic space. The two-sphere $\hat{\mathbb{C}}$ is identified with the
boundary of three-dimensional hyperbolic space $\mathbb{H}^{3}$, and Möbius
transformations are viewed as isometries of $\mathbb{H}^{3}$. The definitions
of crystallographic and Kleinian circle packings are as follows:

######  Definition 2.12.

Let $\mathscr{P}$ denote a collection of circles whose interiors are disjoint
and dense in $\hat{\mathbb{C}}$. $\mathscr{P}$ is a crystallographic packing
if its residual set is the limit set of a geometrically finite reflection
group of isometries of $\mathbb{H}^{3}$. It is Kleinian if its residual set is
the limit set of any geometrically finite group of isometries of
$\mathbb{H}^{3}$.

All polyhedral packings are crystallographic, and all crystallographic
packings are Kleinian.

The following theorem gives sufficient conditions for a packing $\mathscr{P}$
to be crystallographic or Kleinian.

######  Theorem 2.13.

A packing $\mathscr{P}$ satisfying Definition 2.7 is Kleinian if any of the
following conditions hold:

  1. (1)

$B\cup\hat{B}$ is finite.

  2. (2)

After applying some Möbius transformation, $B\cup\hat{B}$ is a strip
configuration (i.e. it contains two parallel lines) and $\Gamma_{2}$ contains
a translation.

  3. (3)

After applying some Möbius transformation, $\Gamma_{2}$ contains two linearly
independent translations.

Furthermore, $\mathscr{P}$ is crystallographic if any of the above conditions
hold, with $\Gamma_{2}$ replaced by the maximal reflective subgroup of
$\Gamma_{2}$.

###### Proof.

By Propositions 2.9, 2.10, the interiors of circles in $\mathscr{P}$ are
disjoint and dense in $\hat{\mathbb{C}}$. We will show that each condition
(1)-(3) implies that $\Gamma$, viewed as a group of isometries of
$\mathbb{H}^{3}$, is geometrically finite. This is essentially a consequence
of Theorems 4.2, 4.10. Up to finite index, we have the following fundamental
domains for $\Gamma_{2}$: in case (1) all of $\hat{\mathbb{C}}$, in case (2) a
strip, and in case (3) a compact parallelogram. Each condition implies that
this fundamental domain only intersects finitely many dual circles in
$\hat{B}$. Using the half-space model of $\mathbb{H}^{3}$, we may define a
half-plane for each wall of the fundamental domain for $\Gamma_{2}$, and a
hemisphere for each dual circle in this fundamental domain. These walls bound
a geometrically finite fundamental domain for $\Gamma$. The limit set of
$\Gamma$ is the residual set of $\mathscr{P}$. Thus the packing is Kleinian.
The proof in the crystallographic case works similarly. ∎

Our construction gives rise to packings which are polyhedral (in fact, all
polyhedral packings), packings which are crystallographic but not polyhedral,
packings which are Kleinian but not crystallographic, and packings which are
none of the above. In Section 6, the wallpaper groups generated by reflections
give rise to crystallographic packings; the wallpaper groups containing
translations but no reflections give rise to Kleinian packings. And the same
refinement method in the proof of Theorem 6.2 can produce aperiodic packings,
which are not Kleinian. These examples can be constructed to have integrality
or superintegrality properties.

Note that Proposition 2.11 and Theorem 2.13 give sufficient but not necessary
conditions for $\mathscr{P}$ to be polyhedral, crystallographic, or Kleinian.
Their proofs involve the structure of the group
$\Gamma=\textsc{Sym}(\mathscr{P},\hat{\mathscr{P}})$. But $\mathscr{P}$ may
have additional symmetry not detected by $\Gamma$; in general,
$\textsc{Sym}(\mathscr{P},\hat{\mathscr{P}})\neq\textsc{Sym}(\mathscr{P})$. We
might make different choices of base and dual configuration $B$, $\hat{B}$,
which give rise to the same packing $\mathscr{P}$ but a different dual packing
$\hat{\mathscr{P}}$, and thus a different $\Gamma$. For example, working with
the classical Apollonian packing, we could select any collection of circles
defining a 3-connected subgraph of the full tangency graph as $B$. With a
nonstandard choice of base and dual configuration, we would find a smaller
symmetry group $\Gamma$, and the packing would not be immediately identifiable
as polyhedral, crystallographic, or Kleinian.

##  3\. Existence and Uniqueness of Packings and a Generalized Koebe-Andreev-
Thurston Theorem

In this section, we recall different versions of the Koebe-Andreev-Thurston
theorem that imply the existence and uniqueness up to conformal automorphism
of our packings in many cases. We conjecture a generalization which would
imply existence and uniqueness in all cases.

To state the problem precisely: let $G$ and $\hat{G}$ be a pair of simple
3-connected plane graphs corresponding to the 1-skeleton of a cellular
decomposition of the sphere $\hat{\mathbb{C}}$ or the plane $\mathbb{C}$, and
the 1-skeleton of the dual cellular decomposition. Note that a cellular
decomposition of the plane can be viewed as a cellular decomposition of the
sphere with a unique accumulation point at $\infty$. Do there exist circle
configurations $B$ and $\hat{B}$, satisfying the conditions of Definition 2.1,
such that $G_{B}\cong G$ and $G_{\hat{B}}\cong\hat{G}$? The $\cong$ symbol
here means an isomorphism of graphs and of the associated cellular
decompositions. Moreover, are the configurations $B$ and $\hat{B}$ unique up
to conformal automorphism of $\hat{\mathbb{C}}$ or $\mathbb{C}$? A conformal
automorphism of $\hat{\mathbb{C}}$ is a Möbius transformation; a conformal
automorphism of $\mathbb{C}$ is a Möbius transformation which fixes infinity,
i.e. a similarity $z\mapsto az+b$ or $z\mapsto a\bar{z}+b$.

The Koebe-Andreev-Thurston theorem answers these questions in the affirmative
when $G$ and $\hat{G}$ are finite graphs. A version of the theorem closely
aligned with this article appears in [5]. Restated in our language, their
Theorem 6 is as follows:

######  Theorem 3.1.

Let $G$ and $\hat{G}$ be a pair of finite, simple, 3-connected plane graphs,
corresponding to a cellular decomposition of $\hat{\mathbb{C}}$ and its dual.
Then there exist circle configurations $B$, $\hat{B}$, satisfying Definition
2.1, such that $G_{B}\cong G$ and $G_{\hat{B}}\cong\hat{G}$, and these
configurations are unique up to Möbius transformation.

As discussed in the introduction, this theorem has been extended to some
infinite graphs. Theorem 4.3, the “Discrete Uniformization Theorem” in [23]
implies the following:

######  Theorem 3.2.

Let $G$ be an infinite simple, 3-connected plane graph, corresponding to a
triangulation of $\mathbb{C}$. Then there exist circle configurations $B$,
$\hat{B}$, satisfying Definition 2.1, such that $G_{B}\cong G$, and these
configurations are unique up to similarity.

Stephenson’s statement of the theorem includes triangulations of other
surfaces as well. The restriction to triangulations, not general cellular
decompositions, seems to be a convenient simplification rather than an
essential restriction. In this case, the existence and uniqueness of $B$
immediately imply the existence and uniqueness of $\hat{B}$. But if $G$ is not
a triangulation, then $B$ and $\hat{B}$ must be constructed together, and both
are needed to ensure uniqueness.

Given these two versions of the Koebe-Andreev-Thurston theorem, it is natural
to conjecture the following common extension:

######  Conjecture 3.3.

Let $G$ and $\hat{G}$ be a pair of simple, 3-connected plane graphs,
corresponding to a cellular decomposition of $\mathbb{C}$ and its dual. Then
there exist circle configurations $B$, $\hat{B}$, satisfying Definition 2.1,
such that $G_{B}\cong G$ and $G_{\hat{B}}\cong\hat{G}$, and these
configurations are unique up to similarity.

This would imply the existence and uniqueness of the whole class of packings
we study. We hope to prove this conjecture in future work, following the
methods of [23]. As further evidence for the conjecture, we remark that
several special cases have been studied carefully in the literature. When $G$
is the triangular lattice, the uniqueness of the associated circle
configuration is a crucial step in Rodin and Sullivan’s celebrated proof of
the convergence of circle packings to the Riemann mapping [17, Appendix 1].
This means that the triangular and hexagonal packings studied in Section 5 are
unique up to similarity. When $G$ is the square lattice, the uniqueness of the
associated circle configuration is the main theorem of [21]. So the square
packing in Section 5 is also unique up to similarity.

It should also be possible to extend Conjecture 3.3 to the hyperbolic plane
and other surfaces, but we have not investigated this.

##  4\. Group Structure Theorems

In this section, we give a complete algebraic description of the symmetry
groups associated to packings $\mathscr{P}$. We begin with the structure of
the groups $\Gamma_{1}$ and $\Gamma_{2}$. Then, via an examination of the
action on $\mathscr{P}$, we show that $\Gamma=\Gamma_{1}\rtimes\Gamma_{2}$.

######  Theorem 4.1.

$\Gamma_{1}$ is a free Coxeter group generated by $\sigma_{d}$ for
$d\in\hat{B}$, where the only relations are $\sigma_{d}^{2}=1$.

###### Proof.

We must show that the relations $\sigma_{d}^{2}=1$ are the only ones. Choose a
point $x\in\hat{\mathbb{C}}$, outside of every circle in $\hat{B}$. Suppose
that a string $\sigma_{d_{1}}\cdots\sigma_{d_{k}}\in\Gamma_{1}$, with
consecutive $d_{i}$ distinct, is applied to $x$. Since $x$ is outside of
$d_{k}$, $\sigma_{d_{k}}(x)$ is in the interior of $d_{k}$. Repeating this
process, by property (1) from Definition 2.1, we find that
$\sigma_{d_{1}}\cdots\sigma_{d_{k}}(x)$ is in the interior of $d_{1}$. Thus it
is not equal to $x$, and $\sigma_{d_{1}}\cdots\sigma_{d_{k}}$ is not the
identity. ∎

######  Theorem 4.2.

If $B$ is finite, then $\Gamma_{2}$ is the group of symmetries of a
polyhedron. If $B$ is infinite, then $\Gamma_{2}$ is conjugate to a discrete
group of isometries of the plane: a cyclic group, dihedral group, frieze
group, or wallpaper group.

A frieze group is a discrete group of isometries of the plane which contains
translations in one direction; a wallpaper group is a discrete group of
isometries of the plane which contains translations in two linearly
independent directions.

###### Proof.

If $B$ is finite, then we have shown that $G_{B}$ is the graph of a
polyhedron. By the Koebe-Andreev-Thurston theorem, the circle configuration
$B$ with this graph is unique up to Möbius transformations. Moreover, any
graph automorphism of $G_{B}$ gives rise to a permutation of the circles in
$G_{B}$, which must be realized by a Möbius transformation in $\Gamma_{2}$.
Conversely, any element of $\Gamma_{2}$ determines a graph automorphism of
$G_{B}$. So $\Gamma_{2}\equiv\textsc{Aut}(G_{B})$. There is a three-
dimensional realization of the polyhedron $G_{B}$, called the canonical
embedding, with all edges tangent to the unit sphere, such that every
automorphism of $G_{B}$ is realized as a rigid motion of this polyhedron
preserving the sphere [25, Thm. 4.13]. Thus $\Gamma_{2}$ is the group of
symmetries of a polyhedron, or a finite group of isometries of the sphere.

If $B$ is infinite, then we may apply a Möbius transformation to place the
unique accumulation point of $B$ at $\infty$. All symmetries of $B$ must map
the accumulation point to itself, so they must have the form $f(z)=az+b$ or
$f(z)=a\bar{z}+b$. If $|a|\neq 1$, then such a map has a fixed point in
$\mathbb{C}$, which will be an attracting fixed point for $f$ or $f^{-1}$.
This produces an additional accumulation point for $B$, a contradiction. Thus
$\Gamma_{2}$ consists of maps $f(z)=az+b$ or $f(z)=a\bar{z}+b$ with $|a|=1$,
which are isometries of the plane $\mathbb{C}$. Again because $B$ has no
accumulation points in $\mathbb{C}$, $\Gamma_{2}$ must be a discrete group of
isometries of $\mathbb{C}$. The rest of the theorem follows from the
classification of these groups, see [6]. ∎

In Section 6, we will show that each of the above possibilities for
$\Gamma_{2}$ is in fact realized by an appropriate choice of $B$.

In order to understand the interactions of $\Gamma_{1}$ and $\Gamma_{2}$, we
need further geometric information about their action on packings. Any circle
in $\mathscr{P}$ is $c=\sigma_{d_{1}}\cdots\sigma_{d_{k}}(c_{0})$ for some
$c_{0}\in B$, $d_{1},\ldots d_{k}\in\hat{B}$. Define the height
$\mathrm{ht}(c)$ as the minimum $k$ for which such an expression exists. A
circle has height $0$ if and only if it is in $B$.

######  Lemma 4.3.

Suppose $c\in\mathscr{P}$ has $\mathrm{ht}(c)>0$. For $d\in\hat{B}$,
$\mathrm{ht}(\sigma_{d}(c))<\mathrm{ht}(c)$ if and only if $c$ is inside $d$.

###### Proof.

Say that $c=\sigma_{d_{1}}\cdots\sigma_{d_{k}}(c_{0})$ and that this
expression is minimal. Then $c_{0}$ is not orthogonal to $d_{k}$, so it is
outside $d_{k}$. Thus $\sigma_{d_{k}}(c_{0})$ is inside $d_{k}$, and hence
outside $d_{k-1}$. Repeating inductively, we conclude that $c$ is inside
$d_{1}$. Since the circles of $\hat{B}$ have disjoint interiors, $d_{1}$ is
the unique circle in $\hat{B}$ whose interior contains $c$. If we reflect
through $d_{1}$, the height of $c$ will be lowered.

On the other hand, if we reflect through any other circle $d\in\hat{B}$,
$\sigma_{d}(c)$ will be inside $d$. Then
$\mathrm{ht}(c)=\mathrm{ht}(\sigma_{d}\sigma_{d}(c))<\mathrm{ht}(\sigma_{d}(c))$,
so the height of $c$ will be raised. ∎

######  Lemma 4.4.

For any $g\in\Gamma$ and $d\in\hat{B}$, $g(B)$ is either completely inside $d$
or completely outside $d$. (Some circles of $g(B)$ are allowed to be
orthogonal to $d$ in either case.)

###### Proof.

For convenience in this proof, we consider a circle orthogonal to $d$ as being
both inside and outside $d$. Since $d$ is a dual circle and $g$ preserves
$\hat{\mathscr{P}}$, we know that $g^{-1}(d)\in\hat{\mathscr{P}}$. By
Proposition 2.8, we know that for some $\tilde{d}\in\hat{B}$, $g^{-1}(d)$ is
either $\tilde{d}$ or inside $\tilde{d}$. All circles in $B$ are completely
outside $\tilde{d}$ by definition, so they are outside $g^{-1}(d)$ as well.
Applying $g$ to both $B$ and $g^{-1}(d)$, we have that all circles in $g(B)$
are completely inside or completely outside $d$. ∎

######  Lemma 4.5.

Let $g$ be a Möbius transformation, $d$ be a circle, and $\sigma_{d}$ be the
reflection across $d$. We have $g\sigma_{d}g^{-1}=\sigma_{g(d)}$.

###### Proof.

The map $g\sigma_{d}g^{-1}\sigma_{g(d)}$ is a holomorphic Möbius
transformation which fixes the circle $g(d)$, so it is the identity. Thus
$g\sigma_{d}g^{-1}=\sigma_{g(d)}$. ∎

######  Proposition 4.6.

For a circle packing $\mathscr{P}$ with base configuration $B$ and dual
configuration $\hat{B}$, we have
$\Gamma=\left\langle\Gamma_{1},\Gamma_{2}\right\rangle$.

###### Proof.

By definition, we have $\Gamma_{1}\leq\Gamma$. For any $g\in\Gamma_{2}$, we
know $g$ preserves $\hat{\mathscr{P}}$ because it preserves $\hat{B}$, and we
also know $g$ preserves $\mathscr{P}$ because it preserves both $B$ and
$\hat{B}$. Hence, we have $\Gamma_{2}\leq\Gamma$ and thus
$\left\langle\Gamma_{1},\Gamma_{2}\right\rangle\subseteq\Gamma$.

We now want to show that $\Gamma\subseteq\langle\Gamma_{1},\Gamma_{2}\rangle$.
By Lemma 4.4, for any $g\in\Gamma$, we know $g(B)$ is either outside all dual
circles or is completely inside some dual circle.

Case 1: If $g(B)$ is outside all dual circles, then $g(B)\subseteq B$. Let
$G_{B}$ be the tangency graph of $B$, and let $G_{g(B)}$ be the subgraph of
$G_{B}$ that is also the tangency graph of $g(B)$. The faces of $G_{g(B)}$ are
either faces of $G_{B}$ or unions of these faces. We will show that each face
of $G_{g(B)}$ is a face of $G_{B}$, implying that $g(B)=B$. Notice that every
face in $G_{g(B)}$ corresponds to a circle $g(\tilde{d})\in\hat{\mathscr{P}}$
for some $\tilde{d}\in\hat{B}$. Since $g(\tilde{d})$ is ringed by circles from
$g(B)$ and thus from $B$, we know $g(\tilde{d})$ cannot be inside any dual
circle. By Proposition 2.8, then, we know that $g(\tilde{d})$ is a dual
circle, which implies that $g(B)=B$. Because $g$ preserves $G_{B}$, we know
$g$ also preserves the dual graph of $G_{B}$, which is the tangency graph of
$g(\hat{B})$. Let $G_{\hat{B}}$ be the dual graph of $G_{B}$. Since there is a
unique circle orthogonal to a ring of circles, every vertex in $G_{\hat{B}}$
must correspond to a circle in $\hat{B}$. Therefore, we have
$g(\hat{B})=\hat{B}$ and thus $g\in\Gamma_{2}$.

Case 2: If $g(B)$ is completely inside some dual circle $d_{1}\in\hat{B}$,
reflect across $d_{1}$ and apply Lemma 4.4 to the configuration
$\sigma_{d_{1}}g(B)$. If this configuration is outside all dual circles, then
by Case 1, we conclude that $\sigma_{d_{1}}g\in\Gamma_{2}$. Otherwise, it lies
inside some other dual circle $d_{2}$, and we can reflect across $d_{2}$ and
repeat the argument. It suffices to show that this process eventually
terminates, i.e. that $\sigma_{d_{k}}\cdots\sigma_{d_{1}}g(B)$ is outside all
the dual circles for some $k$.

Consider a finite set $S$ of circles in $g(B)$ such that no circle in
$\hat{\mathbb{C}}$ is orthogonal to all of them. Such a set can be obtained
starting from three circles ringing a common dual circle, and then choosing a
fourth which is not part of this ring. At each step of reflecting through a
dual circle, the height of each circle in $S$ decreases or stays constant by
Lemma 4.3. Moreover, the latter possibility can only occur if a circle is
orthogonal to the dual, so at least one circle’s height decreases at each
step. The process terminates when all the circles in $S$ reach height $0$, so
they are mapped to $B$. When this occurs, the circles in $S$ are all
orthogonal to or outside each circle in $\hat{B}$, and at least one is outside
each dual circle. By Lemma 4.4, all circles in $g(B)$ are then outside or
orthogonal to each dual circle. ∎

######  Corollary 4.7.

For any $g\in\Gamma$, there exists $g_{1}\in\Gamma_{1}$ and
$g_{2}\in\Gamma_{2}$ such that $g=g_{1}g_{2}$.

###### Proof.

By the proof of Theorem 4.6, every $g\in\Gamma$ is of the form
$g=\sigma_{d_{1}}\cdots\sigma_{d_{k}}g_{2}$ for some sequence of
$d_{i}\in\hat{B}$ and $g_{2}\in\Gamma_{2}$. Let
$g_{1}=\sigma_{d_{1}}\cdots\sigma_{d_{k}}$. By definition, we know that
$g_{1}\in\Gamma_{1}$, and $g=g_{1}g_{2}$. ∎

######  Proposition 4.8.

$\Gamma_{1}$ is a normal subgroup of $\Gamma$.

###### Proof.

By Proposition 4.6, it suffices to show that $\Gamma_{2}$ normalizes
$\Gamma_{1}$. Let $g\in\Gamma_{2}$, and let $\sigma_{d}$ be a generator of
$\Gamma_{1}$. By Lemma 4.5, $g\sigma_{d}g^{-1}=\sigma_{g(d)}$. Since
$g(d)\in\hat{B}$, this is an element of $\Gamma_{1}$. Thus $\Gamma_{2}$, and
hence $\Gamma$, normalizes $\Gamma_{1}$. ∎

In general, $\Gamma_{2}$ is not a normal subgroup of $\Gamma$. Let $g_{2}$ be
an element in $\Gamma_{2}$ that sends some dual circle $d\in\hat{B}$ to a
different dual circle in $\hat{B}$. The element
$\sigma_{d}g_{2}\sigma_{d}^{-1}\in\Gamma$ sends $d$ to a circle inside $d$,
which means that $\sigma_{d}g_{2}\sigma_{d}^{-1}$ doesn’t preserve $\hat{B}$
and thus $\sigma_{d}g_{2}\sigma_{d}^{-1}\not\in\Gamma_{2}$.

######  Proposition 4.9.

The intersection of $\Gamma_{1}$ and $\Gamma_{2}$ is trivial.

###### Proof.

Let $\sigma_{d_{1}}\cdots\sigma_{d_{k}}\in\Gamma_{1}$ with consecutive $d_{i}$
distinct. Choose a dual circle $d\neq d_{k}$ in $\hat{B}$. As in the proof of
Prop. 2.8, the map $\sigma_{d_{1}}\cdots\sigma_{d_{k}}$ sends $d$ to a circle
inside $d_{1}$, so it doesn’t preserve $\hat{B}$ and thus is not in
$\Gamma_{2}$. ∎

As a a direct result of Proposition 4.6, Proposition 4.8, and Proposition 4.9,
we have the following theorem:

######  Theorem 4.10.

$\Gamma\cong\Gamma_{1}\rtimes\Gamma_{2}$.

We conclude this section by sketching some results on the structure of the
supergroup. Let $\hat{\Gamma}_{1}$ denote the group generated by reflections
across the base circles $b\in B$. Recall that the superpacking $\mathscr{S}$
is defined as the orbit of $B$ under the group
$\langle\Gamma_{1},\hat{\Gamma}_{1}\rangle$. The dual superpacking
$\hat{\mathscr{S}}$ is the orbit of $\hat{B}$ under this group. Define
$\Gamma_{\mathscr{S}}=\textsc{Sym}(\mathscr{S},\hat{\mathscr{S}})$, the group
of Möbius transformations which preserve both the superpacking and the dual
superpacking. One has the following:

######  Theorem 4.11.

The group $\langle\Gamma_{1},\hat{\Gamma}_{1}\rangle$ is a Coxeter group with
generators $\sigma_{c}$ for $c\in B$, $\sigma_{d}$ for $d\in\hat{B}$, and
relations $\sigma_{c}^{2}=1$ for all $c\in B$, $\sigma_{d}^{2}=1$ for all
$d\in\hat{B}$, and $\sigma_{c}\sigma_{d}=\sigma_{d}\sigma_{c}$ for all pairs
$c\in B$, $d\in\hat{B}$ intersecting orthogonally.

######  Theorem 4.12.

$\Gamma_{\mathscr{S}}\cong\langle\Gamma_{1},\hat{\Gamma}_{1}\rangle\rtimes\Gamma_{2}$.

We omit the proofs because they are similar to previous ones in this section.

##  5\. Examples

This section introduces examples of our construction. We call our three main
examples the triangular, square, and hexagonal packings. We will focus on the
arithmetic properties of these examples: quadratic forms and linear relations
satisfied by the curvatures, integrality and superintegrality. A “typical”
packing satisfying our definition will have little arithmetic interest, but
the highly symmetric nature of these examples adds more structure.

Some of these examples have appeared in the literature in other contexts, but
their properties have not been explored in detail. The square packing is
discussed in [16, Figure 10.17] as the limit of the $1/n$ cusp groups in
Maskit’s slice. As discussed in the introduction, the triangular and hexagonal
packings appear in Kontorovich and Nakamura’s classification of superintegral
crystallographic packings. And our examples are closely related to the limit
packings in [2]. The limit of pyramid packings is the original Apollonian
packing, the limit of prism packings is our square packing, and the limit of
antiprism packings is our triangular packing. Other families have more
complicated limits, which will be briefly discussed at the end of this
section.

Every packing has linear and quadratic forms satisfied by the curvatures, like
the Descartes quadratic form for the Apollonian packing. The following
definition characterizes packings with number-theoretic structure:

######  Definition 5.1.

A packing $\mathscr{P}$ is integral if every circle in $\mathscr{P}$ has
integral curvature. The packing is superintegral if every circle in the
superpacking has integral curvature.

We will also say that an equivalence class of packings under Möbius
transformations is (super)integral if one packing in the class has this
property.

The main tool to find the linear and quadratic forms, and check
(super)integrality, is an inversive coordinate system for oriented generalized
circles in $\hat{\mathbb{C}}$. A circle is represented as
$(\tilde{b},b,h_{1},h_{2})^{T}\in\mathbb{R}^{4}$, where $b$ is the signed
curvature, $\tilde{b}$ is the curvature after inversion through the unit
circle, and $(h_{1},h_{2})$ are the coordinates of the center, multiplied by
the curvature. Every circle satisfies the quadratic equation
$h_{1}^{2}+h_{2}^{2}-b\tilde{b}=1$. The action of Möbius transformations on
generalized circles becomes a linear action preserving the quadratic form in
this coordinate system. This setup is well explained in [15] and in [11]. The
article [1] gives a full set of linear and quadratic forms for all polyhedral
packings.

In each of the following examples, we begin with base and dual circle
configurations. We check integrality and superintegrality using the inversive
coordinate system. We give quadratic and linear relations sufficient to
determine the curvatures of all the circles in the packing from a finite set
of base circles (in fact, just three). The proofs of these relations are
omitted because they are similar to the polyhedral case. We also relate these
packings to others with commensurate symmetry groups.

###  5.1. Square Packing

The base and dual configurations for this packing are shown in Figure 2 (A).
The dual configuration is a translation of the base configuration. In these
configurations, all the circles can be represented with coordinates
$(\tilde{b},b,h_{1},h_{2})^{T}\in\mathbb{Z}^{4}$, so the packing is
superintegral.

Quadratic and linear relations satisfied by curvatures in the square packing
are shown in Figure 5. In the formulas, $b_{i}$ represents the curvature of
circle $i$. Any image of one of these configurations under $\Gamma$ will
satisfy the same relation. The symmetry group of the square superpacking is
commensurate to the symmetry group of the Apollonian superpacking, as
illustrated in Figure 6. The full packing is shown in Figure 2 (B).

1234$(b_{1}-3b_{2})^{2}+(b_{4}-3b_{3})^{2}=2(b_{1}+b_{2})(b_{3}+b_{4})$

1324$b_{1}+b_{3}=b_{2}+b_{4}$

5678$b_{5}+b_{7}=b_{6}+b_{8}$

Figure 5. Quadratic form and linear relations

(a) The orbit of the blue circle under the group generated by reflections
across all the circles is the square superpacking.

(b) The orbit of the blue circles under the group generated by reflections
across all the circles is the Apollonian superpacking.

Figure 6. Commensurability of superpacking groups

###  5.2. Triangular and Hexagonal Packings

The base configuration and dual configurations for the triangular packing are
shown in Figure 7. In the hexagonal packing, the roles of base and dual
configurations are reversed. Circles in the base and dual configuration for
the triangular packing can be represented with coordinates in the sets

$\displaystyle\left\\{(\tilde{b},b,h_{1},h_{2})^{T}\in\mathbb{R}^{4}\,\left|\,b,\tilde{b}\in\mathbb{Z},\,h_{1}+h_{2}i\in
2\mathbb{Z}\left[\tfrac{1+i\sqrt{3}}{2}\right]\right.\right\\},$
$\displaystyle\left\\{(\tilde{b},b,h_{1},h_{2})^{T}\in\mathbb{R}^{4}\,\left|\,b,\tilde{b}\in\sqrt{3}\mathbb{Z},\,h_{1}+h_{2}i\in
2i\mathbb{Z}\left[\tfrac{1+i\sqrt{3}}{2}\right]\setminus
2\sqrt{3}\mathbb{Z}\left[\tfrac{1+i\sqrt{3}}{2}\right]\right.\right\\}$

respectively. Reflection across the base and dual circles preserve these sets,
so the triangular packing is superintegral, and rescaling by $\sqrt{3}$, we
see that the hexagonal packing is superintegral as well.

Figure 7. Base and dual configurations for the triangular and hexagonal
packings

Quadratic and linear relations satisfied by curvatures in the triangular and
hexagonal packings are shown in Figure 8. The symmetry group of the triangular
and hexagonal superpackings is commensurate to the symmetry group of the limit
of trapezohedral superpackings from [2], as illustrated in Figure 9. We
conjecture that this group is not commensurate to the symmetry group of any
polyhedral superpacking. The full triangular packing is shown in Figure 10 and
the full hexagonal packing is shown in Figure 11. Note: the circles in Figures
2 (B), 10, and 11 were obtained by reflecting a subset of base circles across
four generations of dual circles.

1234$(3b_{1}-b_{2}+3b_{3}-b_{4})^{2}=12b_{1}b_{3}+4b_{2}b_{4}$

1234$2b_{1}^{2}+2b_{2}^{2}+2b_{3}^{2}+10b_{4}^{2}=(b_{1}+b_{2}+b_{3}+b_{4})^{2}$

1234$2b_{1}-2b_{2}=b_{4}-b_{3}$

Figure 8. Quadratic forms and linear relation

(a) The orbit of the blue circle is the triangular superpacking. The orbit of
the red circle is the hexagonal superpacking.

(b) The orbit of the red circles is the limit of trapezohedral superpackings.

Figure 9. Commensurability of superpacking groups

Figure 10. Triangular packing

Figure 11. Hexagonal packing

We finish this section with some brief remarks on the limits of polyhedral
packing families studied in [2]. The limit of trapezohedral packings is
superintegral. As discussed above, the symmetry group of its superpacking is
commensurate to that of the triangular and hexagonal superpackings. The limit
of the cupola packings is also superintegral. The symmetry group of its
superpacking is commensurate to that of the octahedral superpacking. The limit
of the anticupola packings is not integral; this can be proven by an infinite
descent argument.

##  6\. Wallpaper Groups

The goal of this section is to illustrate the rich variety of circle packings
satisfying Definition 2.7, with a focus on symmetry groups. We give the
following converse to Theorem 4.2:

######  Theorem 6.1.

All the possible groups listed in Theorem 4.2 actually arise as the group of
symmetries $\Gamma_{2}=\textsc{Sym}(B,\hat{B})$ for some base and dual
configurations $B$, $\hat{B}$ satisfying Definition 2.1.

In the case of finite $B$, $\hat{B}$, this theorem follows directly from the
Koebe-Andreev-Thurston theorem. For infinite $B$, $\hat{B}$, the wallpaper
group case is Theorem 6.2 below. The other cases are simpler, and the proofs
are omitted.

######  Theorem 6.2.

Any wallpaper group is the symmetry group $\Gamma_{2}$ of the base and dual
configurations $B$, $\hat{B}$ of some circle packing. Moreover, such $B$ can
be realized as the refinement of the base configuration of either the
triangular or square packing.

###### Proof.

We will prove this theorem by illustration. We label the centers of rotations
and axes of reflections using notation introduced in Figure 12. The entire
base configurations can be generated by these symmetries or translations,
starting from the circles shown.

Axis of reflectionAxis of glide reflectionCenter of rotation by $\pi$Center of
rotation by $2\pi/3$Center of rotation by $\pi/2$Center of rotation by $\pi/3$
Figure 12. Labels for axes of reflection, centers of rotation

The wallpaper groups $p6m$ and $p4m$ are the most complicated ones, and they
are the symmetry groups of the base configurations of triangular and square
packings, respectively, as shown in Figure 13.

Refining these two base configurations by adding smaller circles removes
symmetries. Thus, we can obtain base configurations with different symmetry
groups by refining the triangular or square configurations. Notice that the
refinements given below are well-defined base configurations because the
associated dual configurations still exist. Refinements of the triangular
configuration are shown in Figure 14 and refinements of the square
configuration are shown in Figure 15. Some wallpaper groups could be obtained
by refining either configuration, but we only show one realization of each
group.

Wallpaper Group $p6m$ | Wallpaper Group $p4m$
---|---
|
Figure 13. Maximal wallpaper groups

Wallpaper Group $p3$ | Wallpaper Group $p31m$
---|---
|

Wallpaper Group $p3m1$ | Wallpaper Group $p6$
---|---
|

Figure 14. Refinements of the triangular base configuration

Wallpaper Group $p1$ | Wallpaper Group $p2$ | Wallpaper Group $pm$
---|---|---
|  |

Wallpaper Group $pg$ | Wallpaper Group $cm$ | Wallpaper Group $pmm$
---|---|---
|  |

Wallpaper Group $pmg$ | Wallpaper Group $pgg$ | Wallpaper Group $cmm$
---|---|---
|  |

Wallpaper Group $p4$ | Wallpaper Group $p4g$
---|---
|

Figure 15. Refinements of the square base configuration

∎

## References

  * [1] Nooria Ahmed, William Ball, Ellis Buckminster, Emilie Rivkin, Dylan Torrance, Jake Viscusi, Runze Wang, Ian Whitehead, and S. Yang, _Domains of convergence for polyhedral packings_ , arxiv:2109.01289.
  * [2] Nooria Ahmed, William Ball, Emilie Rivkin, Dylan Torrance, Jake Viscusi, and Ian Whitehead, _Fractal dimensions of polyhedral packings_ , unpublished.
  * [3] E. M. Andreev, _Convex polyhedra of finite volume in Lobačevskiĭ space_ , Mat. Sb. (N.S.) 83 (125) (1970), 256–260. MR 0273510
  * [4] Alan F. Beardon and Kenneth Stephenson, _The uniformization theorem for circle packings_ , Indiana Univ. Math. J. 39 (1990), no. 4, 1383–1425. MR 1087197
  * [5] Graham R. Brightwell and Edward R. Scheinerman, _Representations of planar graphs_ , SIAM J. Discrete Math. 6 (1993), no. 2, 214–229. MR 1215229
  * [6] H. S. M. Coxeter, _Introduction to geometry_ , Wiley Classics Library, John Wiley & Sons, Inc., New York, 1989, Reprint of the 1969 edition. MR 990644
  * [7] Elena Fuchs, Katherine E. Stange, and Xin Zhang, _Local-global principles in circle packings_ , Compos. Math. 155 (2019), no. 6, 1118–1170. MR 3952499
  * [8] Ronald L. Graham, Jeffrey C. Lagarias, Colin L. Mallows, Allan R. Wilks, and Catherine H. Yan, _Apollonian circle packings: number theory_ , J. Number Theory 100 (2003), no. 1, 1–45. MR 1971245
  * [9] Gerhard Guettler and Colin Mallows, _A generalization of Apollonian packing of circles_ , J. Comb. 1 (2010), no. 1, [ISSN 1097-959X on cover], 1–27. MR 2675919
  * [10] Michael Kapovich and Alex Kontorovich, _On superintegral kleinian sphere packings, bugs, and arithmetic groups_ , arxiv:2104.13838.
  * [11] Jerzy Kocik, _A theorem on circle configuarations_ , arXiv:0706.0372.
  * [12] Paul Koebe, _Kontaktprobleme der konformen Abbildung_ , Ber. Verh. Sächs. Akad. Leipzig 88 (1936), 141–164.
  * [13] Alex Kontorovich and Kei Nakamura, _Geometry and arithmetic of crystallographic sphere packings_ , Proc. Natl. Acad. Sci. USA 116 (2019), no. 2, 436–441. MR 3904690
  * [14] Alex Kontorovich and Hee Oh, _Apollonian circle packings and closed horospheres on hyperbolic 3-manifolds_ , J. Amer. Math. Soc. 24 (2011), no. 3, 603–648, With an appendix by Oh and Nimish Shah. MR 2784325
  * [15] Jeffrey C. Lagarias, Colin L. Mallows, and Allan R. Wilks, _Beyond the Descartes circle theorem_ , Amer. Math. Monthly 109 (2002), no. 4, 338–361. MR 1903421
  * [16] David Mumford, Caroline Series, and David Wright, _Indra’s pearls: the vision of Felix Klein_ , Cambridge University Press, Cambridge, 2015. MR 3558870
  * [17] Burt Rodin and Dennis Sullivan, _The convergence of circle packings to the Riemann mapping_ , J. Differential Geom. 26 (1987), no. 2, 349–360. MR 906396
  * [18] Rudolf Scharlau and Claudia Walhorn, _Integral lattices and hyperbolic reflection groups_ , no. 209, 1992, Journées Arithmétiques, 1991 (Geneva), pp. 15–16, 279–291. MR 1211022
  * [19] Oded Schramm, _Existence and uniqueness of packings with specified combinatorics_ , Israel J. Math. 73 (1991), no. 3, 321–341. MR 1135221
  * [20] by same author, _Rigidity of infinite (circle) packings_ , J. Amer. Math. Soc. 4 (1991), no. 1, 127–149. MR 1076089
  * [21] by same author, _Circle patterns with the combinatorics of the square grid_ , Duke Math. J. 86 (1997), no. 2, 347–389. MR 1430437
  * [22] Katherine E. Stange, _The Apollonian structure of Bianchi groups_ , Trans. Amer. Math. Soc. 370 (2018), no. 9, 6169–6219. MR 3814328
  * [23] Kenneth Stephenson, _Introduction to circle packing_ , Cambridge University Press, Cambridge, 2005, The theory of discrete analytic functions. MR 2131318
  * [24] W.P. Thurston, _The geometry and topology of three-manifolds_ , Princeton University, 1979.
  * [25] Günter M. Ziegler, _Lectures on polytopes_ , Graduate Texts in Mathematics, vol. 152, Springer-Verlag, New York, 1995. MR 1311028


# A Parallel SystemC Virtual Platform for Neuromorphic Architectures

Melvin Galicia, Farhad Merchant and Rainer Leupers  Institute for
Communication Technologies and Embedded Systems
RWTH Aachen University
Aachen, Germany
{galicia, farhad.merchant<EMAIL_ADDRESS>
###### Abstract

With the increasing interest in neuromorphic computing, designers of embedded
systems face the challenge of efficiently simulating such platforms to enable
architecture design exploration early in the development cycle. Executing
artificial neural network applications on neuromorphic systems which are being
simulated on virtual platforms (VPs) is an extremely demanding computational
task. Nevertheless, it is a vital benchmarking task for comparing different
possible architectures. Therefore, exploiting the multicore capabilities of
the VP’s host system is essential to achieve faster simulations. Hence, this
paper presents a parallel SystemC based VP for RISC-V multicore platforms
integrating multiple computing-in-memory neuromorphic accelerators. In this
paper, different VP segmentation architectures are explored for the
integration of neuromorphic accelerators and are shown their corresponding
speedup simulations compared to conventional sequential SystemC execution.

###### Index Terms:

Electronic System Level, Parallel SystemC Simulation, Neuromorphic Systems,
Computing-in-Memory

##  I Introduction

The limitations of the von Neumann model have forced researchers to look into
alternate computer architectures. The von Neumann model fares poorly due to
the energy spent in data movements especially for machine learning
applications in embedded edge devices, where a huge amount of unstructured
data is to be processed [1]. Emerging non-volatile memory technology such as
resistive random access memory (ReRAM) offers a unique advantage in-terms of
computing-in-memory features along with storage capabilities [2]. Such a
feature alleviates the need for data movement between a processor and memory
subsystem significantly, resulting in energy efficiency and high-performance
of a neuromorphic architecture. The design and development of a neuromorphic
platform requires careful consideration of several architectural parameters
[3]. In addition, the software’s performance to be executed needs to be tuned
to attain the best possible performance at an early stage in the design cycle.

A virtual platform (VP) offers the simulation of architectural components to
achieve the right balance between time to market and the performance of the
target product [4]. Especially for complex platforms, VPs offer a unique
advantage that the software development can begin significantly earlier than
the hardware prototype’s availability [5].

The emergence of neuromorphic computing-in-memory (CIM) architectures in the
last decade has propelled research to develop various simulation platforms for
performance prediction and early-stage software development. While some of
these simulators lack the required accuracy, others are extremely slow due to
the details involved in the simulation. Additionally, the application
benchmarks used in machine learning are extremely demanding computational
tasks, which mainly require vector-matrix multiplication (VMM) operations that
also have to be executed within the VP. All this creates the need of a VP
capable of effectively using all computing power available to its host system.
So, it can speed up the simulations required to perform investigation of
neuromorphic architectures. To address these issues related to neuromorphic
VPs, this work introduces a SystemC based VP that utilizes, since its
conception, parallelization techniques to efficiently exploit all available
CPU resources in the host system. The major contributions of this paper are as
follows:

  * •

A complete full system-level neuromorphic VP able to run a heavy load of VMM
benchmarks in a parallel simulation environment.

  * •

Simulation segmentation strategies to distribute simulation load among
multiple cores of the host system in a comprehensive method.

  * •

Utilization of time-decoupled parallelization techniques to further enhance
the VP’s simulation speedup.

This paper is structured as follows: Section II provides related work and
background. Section III introduces the structure of the VP and its
characteristics. Section IV summarizes the used time-decoupled and parallel
SystemC methodologies. In section V, the significance of time-decoupled
parallel SystemC VP is demonstrated by using it to perform two different VP
segmentation strategies of multicore and multi-CIM systems. It also presents
benchmark results versus sequential simulations for all analyzed cases.
Finally, in section VI conclusions are drawn and an outlook of future research
is provided.

##  II Related Work

Figure 1: Single-core System-level virtual platform. (a) VP high level architecture, (b) Neuromorphic computing-in-memory architecture: CIM-Unit, (c) Memristor-based crossbar array. TABLE I: Qualitative comparison of the simulators in the literature Simulator |  Architecture-level |  System-level |  Circuit-level |  Exploration possibility |  Parallelization |  CIM Support |  Accelerator-enabled |  Time decoupling-enabled
---|---|---|---|---|---|---|---|---
NVSim [6] |  ✗ |  ✗ |  ✓ |  ✓ |  ✗ |  ✗ |  ✗ |  ✗
NVMain [7] |  ✓ |  ✗ |  ✗ |  ✓ |  ✗ |  ✗ |  ✗ |  ✗
MNSim [8] |  ✓ |  ✗ |  ✗ |  ✗ |  ✗ |  ✓ |  ✗ |  ✗
CorssSim [9] |  ✗ |  ✗ |  ✓ |  ✓ |  ✗ |  ✓ |  ✗ |  ✗
NeuroSim [10] |  ✗ |  ✗ |  ✓ |  ✓ |  ✗ |  ✓ |  ✗ |  ✗
NVM-SPICE [11] |  ✗ |  ✗ |  ✓ |  ✓ |  ✗ |  ✗ |  ✗ |  ✗
AIHWKIT [12] |  ✗ |  ✗ |  ✓ |  ✓ |  ✗ |  ✓ |  ✗ |  ✗
CIMSIM [13] |  ✗ |  ✓ |  ✓ |  ✗ |  ✓ |  ✓ |  ✗ |  ✗
Lee et al. [14] |  ✗ |  ✓ |  ✗ |  ✓ |  ✗ |  ✓ |  ✗ |  ✗
MultiPULPly [15] |  ✗ |  ✓ |  ✗ |  ✓ |  ✓ |  ✓ |  ✗ |  ✗
Proposed work |  ✓ |  ✓ |  ✓ |  ✓ |  ✓ |  ✓ |  ✓ |  ✓

With the advancements in emerging non-volatile memory technologies, various
research groups have tried to explore the technologies at the system level. On
the other hand, for the CMOS-based electronic systems, SystemC-based virtual
prototypes are extremely popular for early software development and
performance tuning. In this section, we first review various simulators and
then briefly discuss the relevant systemC-based virtual prototypes.

A system-level simulation of a tiled neuromorphic architecture is presented in
[13]. The authors have implemented a micro-instruction set to perform the
operations on analog/digital memristive crossbars. Along with the crossbars,
each tile consists of analog/digital converters, digital modulators, and
sample-hold mechanisms. The micro-instruction has three types of functions:
initialize, compute, and read. The system-level simulator presented in [14] is
built on similar principles. NVMain [7] and NVMain 2.0 [16] are the simulators
that can simulate memory and memory interfaces. NVMain focuses on memory-
oriented simulations and requires significant engineering efforts in modeling
and simulation of CIM.

The MNSIM [8] and MNSIM 2.0 [17] simulators use behavioral models to estimate
the worst case and accuracy. The increased simulation speed comes at the cost
of accuracy in MNSIM. Apart from the system level and architecture simulators,
there are various simulators capable of performing circuit-level simulations
[6] [10]. While these simulators can perform low-level simulations and provide
better estimates of the energy efficiency of the design, scaling to system-
level is not possible due to simulation complexity and run-time.

VPs commonly use SystemC as their underlying simulation engine, which deploys
a standard discrete event approach to simulate the concurrency of
multiprocessor systems. There are several approaches for parallelization of
SystemC [18, 19, 20]. A synchronous approach to parallel SystemC simulation is
presented in [20]. Synchronous techniques, however, are unable to exploit the
full parallelism of a simulation, given that they enforce strict time
synchronization among all threads. This synchronous limitation can be overcome
by applying a time-decoupled technique as proposed in [21]. These techniques
have been applied in [22] to multicore VPs for specific workloads like booting
operating systems and for network intensive communication applications.
However for heavy VMM operations needed in machine learning applications, that
are executed in a system-level VP with multiple neuromorphic accelerators
architecture, there is limited work available. Extensive work on different
simulators for neuromorphic architectures with different abstraction levels is
summarized in [1], where the majority of simulators focus on circuit level.
Only a few works like [23, 15] have presented complete integration of
neuromorphic accelerators into a full system. However, they have not been
explicitly developed to exploit the multicore nature of their host systems.
Therefore, this work presents a time-decoupled parallel simulation framework
that allows different VP segmentation strategies of neuromorphic
architectures. A summary for several neuromorphic simulation frameworks is
captured in Table I where a qualitative comparison can be observed. The
proposed simulation framework has a few features that are desired but absent
in simulators found in the literature.

##  III Virtual Platform

The proposed system-level VP consists of several SystemC modules: a set of
multicore processors, L1 instruction and data caches, communication buses, a
main memory module and multiple memristor-based CIM neuromorphic accelerators.
Each SystemC module can work as a transaction initiator and a transaction
target depending on its current function or just as a forwarder in the case of
the bus modules. The modules have their corresponding sockets for TLM2.0
transaction-based communication where a tracing monitoring of each transaction
is being performed and stored for post-processing performance analysis of the
whole system. Fig. 1a depicts the VP high level architecture and illustrates
the tracing approach. By collecting all the transactions, it is possible to
have a histogram and the total number of transactions that occurred during the
execution of any particular benchmark.

###  III-A RISC-V

The VP uses RISC-V [24], an open-source instruction set architecture, as the
main processor of the multicore system. The RISC-V core implemented in [25]
using SystemC and TLM-2.0 is the one integrated into the VP. This
implementation offers a 32/64-bit RISC-V core supporting the IMAC instruction
set. By adding this SystemC RISC-V 64IMAC, the VP achieves a significantly
faster simulation compared to using an RTL implementation. The RISC-V
processor is used to execute different software applications that exploit the
available neuromorphic accelerator by offloading artificial neural network
(ANN) main operations, i.e. VMMs, meanwhile being free to perform other
computational tasks when necessary. It is important to point out that the
RISC-V core does not make use of any standard extension instructions to
offload operations to the neuromorphic accelerator. It is the CIM unit that
interprets and handles the executions of the operations using its internal
controller and own set of micro-instructions as explained in next section.

###  III-B Memristor-Based CIM-Unit

The VP integrates multiple CIM-Units of the one proposed in [13], since it
provides built-in interfaces that allow interaction with SystemC based VPs in
several abstraction levels. Each CIM-Unit consists of two main parts: A
calculator with its digital surrounding and a micro-engine as shown in Fig.
1b. The calculator is a simulator that not only replicates the functional
behavior of a memristor crossbar shown in Fig. 1c, but also covers the
operation of surrounding mixed-signal circuitries, e.g., analog to digital
converter (ADC), digital to analog converter (DAC) and sample and hold (S+H)
elements, which are essential for driving the memristor crossbar. The micro-
engine comprises the digital components, e.g., controller, registers files,
and buffers, that are necessary for operating the memristor crossbar. In
addition to these, the CIM-Unit offers a micro-instruction set that allows the
unit to communicate with the main processor or other functional units, in a
coarse grain reconfigurable architecture. Thus, it can be deployed as an
accelerator as it is the case in this work. For the CIM-Unit to operate it
needs first be supplied with configuration parameters, e.g., size of the
matrix that is to be mapped on the crossbar (a fixed size of 256x256 was
configured), I/O resolution as well as others necessary parameters detailed in
[13]. The controller stores all these parameters which then are passed to the
CIM-Unit through instructions via the configuration register. As soon as all
the parameters are available, the controller, which is a state-machine,
transitions into the state IN signaling that is ready to accept input data. It
stays in IN until all the input data is received. The number of cycles that it
takes to be ready to move to the next state depends on the width of the input
ports, the resolution of the input data, and the size of the input vector. In
the next state, OP, the controller executes the operation that is specified in
the configuration register. In the end, the controller moves to the final
state, OUT, where it sends out the processed information and goes back to the
initial state, IDLE, as soon as all the output is sent out.

###  III-C Memory Modules

The main memory (DRAM), which works as a shared module for all segments, and
instruction and data cache (SRAM) modules, from the work presented in [26],
are the ones integrated in the VP. Since in addition to commonly modeling
memory in SystemC with a read and write access delay, further delays for page
switches and write-to-read switches are introduced.

##  IV Parallel SystemC Simulation

This work implements time-decoupled parallel simulation using techniques
investigated in [27], and also presents different segmentation methodologies.
The segmentation process can be used to explore effective exploitation
configurations for parallel execution of VMM benchmarks in a multicore system
integrating multiple neuromorphic accelerators. Aforementioned implementation
and methodologies are briefly summarized as follows:

Figure 2: SystemC parallel simulation overview.

Figure 3: SytstemC sequential and time-decoupled parallel simulation quantum.

###  IV-A Parallel Simulation Architecture

A complete simulation consists of n segments, defined during the segmentation
process, that are interconnected via dedicated communications channels. By
organizing the simulation as a list of segments each containing its own
SystemC kernel, parallel simulation becomes simple, since it is not necessary
to modify the internal discrete event simulation of each kernel. As a result,
each segment can be stepped independently and is then parallelized by the host
system. However, communication, time-decoupling and post-synchronization
become key procedures during the whole simulation execution. Therefore, the
simulation controller is used to advance the simulation by stepping each
segment and enforcing the sync mechanism. The simulation flow that executes
SystemC for the simulation segments in parallel is shown in Fig. 2.

###  IV-B Time-Decoupled Simulation

Higher simulation speed can be achieved by deploying time-decoupling between
the segments [21]. This is achieved by relaxing the synchronization
requirement between segments and allowing the simulation controller to let the
simulation run in parallel for a certain time before synchronization is
enforced. The controller organises instructions into quanta: a chunk of N
instructions that are grouped together and executed, before simulation time is
advanced by a corresponding amount. The threads runtime are constructed based
on three communication primitives: setup: loads the thread and sets quantum
length to N, exec: asynchronously executes the next N instructions, sync:
blocks until the thread has finished its current quantum. This is illustrated
in Fig. 3 for both cases, when parallel execution is performed and when it is
done sequentially by default in SystemC.

Since each segment now operates in its own time zone with a local time, inter-
segment communication is assigned a latency. Channel latency describes the
amount of time allowed to pass after putting a message into a channel before
it must be fetched by the receiver. It is used by the simulation controller to
determine the amount of time a segment may simulate ahead of time before
risking missing a channel message from a peer segment. Hence, it is the task
of the simulation controller to make sure that no segment has advanced too far
ahead in time. Segments that have not reached their limit time are considered
in simulation, while others are considered waiting for their peers to match
time.

###  IV-C Segmentation

The segmentation process allows the designer to explore different partition
strategies of the VP architecture for the integration of neuromorphic CIM-
Units within a multicore system. There are several parameters that can be
considered to build a complete system, e.g. the total number of CPUs, the
total number of CIM-Units and also the division of load execution. In Fig.
4(a) the segmentation corresponds to an equal number of CIM-Units accessible
to each CPU, i.e. uniform VP segmentation. It is as well possible to segment
in a way that one CPU is dedicated to manage all CIM-Units while other CPUs
are free to perform other tasks, i.e. load-oriented VP segmentation as shown
in Fig. 4(b). Subsequently, the CIM-Units can be allocated in parallel
segments to speed up even further the simulation. Segments form virtual
sequential environments which are allowed to share data using communication
channels. Segments are implemented as shared object files with unique names.
Hence, the system can easily be extended to a n-core system by instantiating
more segments and creating new channel connections. Once the simulation
starts, the controller selects the first segment from the ready to simulate
queue and starts its execution each in its corresponding thread.

##  V Experimental Evaluation

To demonstrate the efficacy of the proposed approach, it is evaluated in
various experiments. The evaluation considers VMM loads from well known ANN
layers to assess the performance benefits of the parallel simulation approach.
Runtimes for all workloads are analyzed using the conventional sequential (sq)
and the parallel simulation approach (pll).

###  V-A Multicore and Multi-CIM Systems

The designs of the VP with uniform and load-oriented segmentation, shown in
Fig. 4(a) and Fig. 4(b) respectively, are explored in this work. In the first
case only two segments are initialized, each is allocated with one RISC-V
processor and with only two CIM-Units. In the second case, four segments are
initialized, segment 0 contains a RISC-V processor and the main shared memory.
Segment 1 contains only one RISC-V processor, meanwhile segment 2 and 3 each
contains only two CIM-Units. The VP components and their parameters are listed
in Table II as well as for the host machine.

###  V-B Benchmarks

Convolutional layers from Googlenet [28], ImageNet [29], and MobileNets [30]
ANNs are selected as benchmarks that run on the VP, these are heavy workloads
that require acceleration for the vector matrix multiplication operations. Due
to the multicore nature of the simulated VP, it is possible to run the
convolutional layers using one CPU and main memory as well as running them by
offloading the operation to the available CIM-Units. The layers are listed in
Table III, where h and w are the matrix height and width respectively and p is
the number of vectors. The matrix and vector sizes varied to assess different
scenarios. Nonetheless, all matrices fit into the configured crossbar inside
the CIM-Units. A simple nested loop algorithm for a VMM operations was
implemented to obtain the O matrix which is the dot product between matrices A
and B. The algorithm was used for the case of RISC-V plus the shared main
memory and also for the case when the host system offloads the operations to
the CIM-Units. However for the latest case, the internal operation handling is
left to the micro-engine component of the CIM-Units. The nested loop algorithm
is described below:

(a)

(b)

(c)

(d)

Figure 4: (a) VP with uniform segmentation, (b) VP with load-oriented
segmentation, (c) VMM benchmarks results for VP with uniform segmentation, and
(d) VMM benchmarks results for VP with load-oriented segmentation.

⬇

1 // Simple VMM algorithm: $O=A(h\times w)\cdot B(w\times p)$

2 for ( i = 1 to h ) {

3 for ( j = 1 to p ) {

4 $O_{(}i,j)=0$

5 for ( k = 1 to w ) {

6 $O_{(}i,j)=O_{(}i,j)+A_{(}i,k)\cdot B_{(}k,j)$

7 } } }

TABLE II: Multicore VP and Host System Components VP Component | Parameters
---|---
RISC-V processor | 2 core 64bits; 1.7GHZ; in order
L1 I and D caches | SRAM 16KB and 32KB
Main memory RAM | DRAM 128MB
CIM-Unit | 2 x Segment ; Crossbar size 256x256
Host Component | Parameters
AMD Ryzen9 processor | 12 core 64bits; 2.2GHZ
L1 I and D caches | 512KB
Main memory RAM | 61GB

###  V-C Performance Results

Results of the simulation runtime of the complete VP systems corresponding to
the segmentation illustrated in Fig. 4(a) and Fig. 4(b) are shown in Fig. 4(c)
and in Fig. 4(d) respectively. In both cases the VP is executing the
benchmarks listed in Table III, once using the conventional sequential (sq)
approach and once the parallel simulation (pll) approach. Additionally, for
the time-decouple technique the quanta of N instructions was set to 10K after
some experimentation. Since increasing the quantum size results in a
significant speedup but not uniformly for all benchmarks. When the quantum
exceeds the above mentioned value the speed decreases due to the higher number
of synchronizations required by the execution of the VMM operations in RISC-V
plus main shared memory. When a sync is needed all instructions inside the
quantum will remain in standby until that sync is performed and consequently
decreasing the achieved speedup. This effect is not seen when the VMM
operations are executed by the CIM-Units due to their extremely low
requirement to access new values inside the main memory. Hence, confirming the
fact that neuromorphic computing, in this case: computing-in-memory,
alleviates the von Neumann bottleneck.

TABLE III: Network Layer Benchmarks Network Name | Layer Type id | h | w | p
---|---|---|---|---
Googlenet | Conv 1 | 224 | 224 | 7
Googlenet | Conv 2 | 56 | 56 | 3
ImageNet | Conv 1 | 224 | 224 | 11
ImageNet | Conv 2 | 207 | 207 | 5
MobileNets | Conv 1 | 224 | 224 | 3
MobileNets | Conv 2 | 112 | 112 | 3

From all cases analyzed for the uniform segmentation an up to 2.3× speedup
efficiency was registered for ImageNet-conv1 and Googlenet-conv1 network
layers. A fast analysis shows that, the total sequential runtime is the sum of
segment 0 and 1 individual runtimes as ilustrated in Fig. 3. Therefore in this
case, it is visible that the obtained speedup is directly related to the
number of system host threads performing the simulation, i.e. two parallel
segments yield at least 2× speedup, when compared to sequential, and the rest
speedup is attributed to the utilization of the time-decoupled technique.
Additionally, from all cases for the load-oriented segmentation an up to 3.3×
speedup efficiency was registered for ImageNet-conv1 and Mobilenets-conv1
network layers. However in this case, the total speedup relation is not so
evident. Again, the total sequential runtime will be the sum of segment 0 to 3
runtimes. Here, from a close examination at the performed load-oriented
segmentation, it is posssible to see that segment 1 does not have a heavy load
and therefore its contribution to the total speedup is minimal. Nonetheless,
the remaining three segments redistribute the total load and yield more than
3× speedup in all benchmarks.

##  VI Conclusion

This work presented a parallel SystemC virtual platform that enables fast
simulations for the integration of neuromorphic accelerators within a
multicore system. Results show that the utilized simulation approach improves
simulation speed by up to 2.3× when a uniform segmentation is performed and up
to 3.3× when it is load-oriented for analyzed workloads in both cases
respectively. Since current host computers offer multicore computational power
to drive several number of threads, exploiting parallel SystemC simulation as
a mean to carry out investigation of VPs for neuromorphic systems is shown to
be an effective methodology.

Future work includes the development of automatic segmentation mechanism when
the load specifications for the neuromorphic accelerators are known
beforehand. Hence, it will allow for a quick starting point to the designer.

## Acknowledgement

This work was supported by BMBF in Germany within NEUROTEC, project numbers
16ES1134 and 16ES1133K.

## References

  * [1] F. Staudigl, F. Merchant, and R. Leupers, “A Survey of Neuromorphic Computing-in-Memory: Architectures, Simulators and Security,” IEEE Design & Test, 2021.
  * [2] A. Lu, X. Peng, and S. Yu, “Compute-in-RRAM with Limited On-chip Resources,” in 2021 IEEE 3rd AICAS, pp. 1–4, 2021.
  * [3] H. H. Li, Y. Chen, C. Liu, J. P. Strachan, and N. Davila, “Looking Ahead for Resistive Memory Technology: A broad perspective on ReRAM technology for future storage and computing,” IEEE Consumer Electronics Magazine, vol. 6, no. 1, pp. 94–103, 2017.
  * [4] O. Bringmann, W. Ecker, A. Gerstlauer, A. Goyal, D. Mueller-Gritschneder, P. Sasidharan, and S. Singh, “The next generation of virtual prototyping: Ultra-fast yet accurate simulation of hw/sw systems,” in DATE, DATE ’15, (San Jose, CA, USA), p. 1698–1707, EDA Consortium, 2015.
  * [5] T. Schmidt, K. Grüttner, R. Dömer, and A. Rettberg, “A program state machine based virtual processing model in systemc,” SIGBED Rev., vol. 11, p. 7–12, Jan. 2015.
  * [6] X. Dong, C. Xu, Y. Xie, and N. P. Jouppi, “Nvsim: A circuit-level performance, energy, and area model for emerging nonvolatile memory,” IEEE TCAD, vol. 31, no. 7, pp. 994–1007, 2012.
  * [7] M. Poremba and Y. Xie, “Nvmain: An architectural-level main memory simulator for emerging non-volatile memories,” in IEEE ISVLSI 2012, pp. 392–397, 2012.
  * [8] L. Xia, B. Li, T. Tang, P. Gu, X. Yin, W. Huangfu, P.-Y. Chen, S. Yu, Y. Cao, Y. Wang, Y. Xie, and H. Yang, “Mnsim: Simulation platform for memristor-based neuromorphic computing system,” in DATE 2016, pp. 469–474, 2016.
  * [9] S. Agarwal, R. B. Jacobs Gedrim, A. H. Hsia, D. R. Hughart, E. J. Fuller, A. A. Talin, C. D. James, S. J. Plimpton, and M. J. Marinella, “Achieving ideal accuracies in analog neuromorphic computing using periodic carry,” in 2017 Symposium on VLSI Technology, pp. T174–T175, 2017.
  * [10] P.-Y. Chen, X. Peng, and S. Yu, “Neurosim: A circuit-level macro model for benchmarking neuro-inspired architectures in online learning,” IEEE TCAD, vol. 37, no. 12, pp. 3067–3080, 2018.
  * [11] H. Yu and Y. Wang, “Nonvolatile State Identification and NVM SPICE,” in Design Exploration of Emerging Nano-scale Non-volatile Memory (H. Yu and Y. Wang, eds.), pp. 45–83, New York, NY: Springer, 2014\.
  * [12] M. J. Rasch, D. Moreda, T. Gokmen, M. L. Gallo, F. Carta, C. Goldberg, K. E. Maghraoui, A. Sebastian, and V. Narayanan, “A flexible and fast pytorch toolkit for simulating training and inference on analog crossbar arrays,” Apr. 2021.
  * [13] A. BanaGozar, K. Vadivel, J. Multanen, P. Jääskeläinen, S. Stuijk, and H. Corporaal, “System Simulation of Memristor Based Computation in Memory Platforms,” in Embedded Computer Systems: Architectures, Modeling, and Simulation, pp. 152–168, Springer International Publishing, 2020\.
  * [14] M. K. F. Lee, Y. Cui, T. Somu, T. Luo, J. Zhou, W. T. Tang, W.-F. Wong, and R. S. M. Goh, “A system-level simulator for rram-based neuromorphic computing chips,” ACM Trans. Archit. Code Optim., vol. 15, Jan. 2019.
  * [15] A. Eliahu, R. Ronen, P.-E. Gaillardon, and S. Kvatinsky, “MultiPULPly: A Multiplication Engine for Accelerating Neural Networks on Ultra-Low-Power Architectures,” J. Emerg. Technol. Comput. Syst., vol. 17, Apr. 2021.
  * [16] M. Poremba, T. Zhang, and Y. Xie, “Nvmain 2.0: A user-friendly memory simulator to model (non-)volatile memory systems,” IEEE Computer Architecture Letters, vol. 14, no. 2, pp. 140–143, 2015.
  * [17] Z. Zhu, H. Sun, K. Qiu, L. Xia, G. Krishnan, G. Dai, D. Niu, X. Chen, X. S. Hu, Y. Cao, Y. Xie, Y. Wang, and H. Yang, MNSIM 2.0: A Behavior-Level Modeling Tool for Memristor-Based Neuromorphic Computing Systems, p. 83–88.  New York, NY, USA: Association for Computing Machinery, 2020.
  * [18] N. Ventroux, J. Peeters, T. Sassolas, and J. C. Hoe, “Highly-parallel special-purpose multicore architecture for SystemC/TLM simulations,” in 2014 International Conference on Embedded Computer Systems: Architectures, Modeling, and Simulation (SAMOS XIV), pp. 250–257, 2014.
  * [19] W. Chen, X. Han, and R. Dömer, “Out-of-order parallel simulation for ESL design,” in DATE 2012, pp. 141–146, 2012.
  * [20] C. Schumacher, J. H. Weinstock, R. Leupers, G. Ascheid, L. Tosoratto, A. Lonardo, D. Petras, and T. Grötker, “legaSCi: Legacy SystemC Model Integration into Parallel Systemc Simulators,” in 2013 IEEE International Symposium on Parallel Distributed Processing, Workshops and Phd Forum, pp. 2188–2193, 2013.
  * [21] J. H. Weinstock, C. Schumacher, R. Leupers, G. Ascheid, and L. Tosoratto, “Time-decoupled parallel SystemC simulation,” in DATE 2014, pp. 1–4, 2014.
  * [22] J. H. Weinstock, R. L. Bücs, F. Walbroel, R. Leupers, and G. Ascheid, “Work-in-Progress: AMVP - A High Performance Virtual Platform using Parallel SystemC for Multicore ARM Architectures,” in 2018 International Conference on Hardware/Software Codesign and System Synthesis (CODES+ISSS), pp. 1–2, 2018.
  * [23] M. Galicia, A. BanaGozar, K. Sturm, F. Staudigl, S. Stuijk, H. Corporaal, and R. Leupers, “NeuroVP: A System-Level Virtual Platform for Integration of Neuromorphic Accelerators,” in IEEE International System-on-Chip Conference,to be published, IEEE SOCC Conference Proceedings, 2021.
  * [24] A. Waterman, Y. Lee, D. A. Patterson, and K. Asanović, “The RISC-V Instruction Set Manual, Volume I: User-Level ISA, Version 2.0,” Tech. Rep. UCB/EECS-2014-54, EECS Department, University of California, Berkeley, May 2014\.
  * [25] V. Herdt, D. Große, P. Pieper, and R. Drechsler, “RISC-V based virtual prototype: An extensible and configurable platform for the system-level,” Journal of Systems Architecture, vol. 109, p. 101756, 2020.
  * [26] G. Onnebrink, R. Leupers, and G. Ascheid, “ESL Black Box Power Estimation: Automatic Calibration for IEEE UPF 3.0 Power Models,” in RAPIDO 2018, RAPIDO ’18, (New York, NY, USA), Association for Computing Machinery, 2018.
  * [27] J. H. Weinstock, R. Leupers, G. Ascheid, D. Petras, and A. Hoffmann, “SystemC-link: Parallel SystemC simulation using time-decoupled segments,” in DATE 2016, pp. 493–498, 2016.
  * [28] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” in 2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 1–9, 2015.
  * [29] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “ImageNet Classification with Deep Convolutional Neural Networks,” Commun. ACM, vol. 60, p. 84–90, May 2017.
  * [30] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam, “MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications,” CoRR, vol. abs/1704.04861, 2017\.