# Sandpile toppling on Penrose tilings:
identity and isotropic dynamics

Jérémy Fersula  Université publique  Camille Noûs  Cogitamus laboratory  Kévin
Perrot  Université publique

###### Abstract

We present experiments of sandpiles on grids (square, triangular, hexagonal)
and Penrose tilings. The challenging part is to program such simulator; and
our javacript code is available online, ready to play! We first present some
identity elements of the sandpile group on these aperiodic structures, and
then study the stabilization of the maximum stable configuration plus the
identity, which lets a surprising circular shape appear. Roundness
measurements reveal that the shapes are not approaching perfect circles,
though they are close to be. We compare numerically this almost isotropic
dynamical phenomenon on various tilings.

## Preamble

The experiments presented in this paper were conducted using JS-Sandpile, a
javascript sandpile simulator we developed. The code is available at

https://github.com/huacayacauh/JS-Sandpile,

and it is ready for the reader to play at

https://huacayacauh.github.io/JS-Sandpile/.

The color codes of all our pictures are a gradation, from almost white = 0
grain, to almost black = $deg(v)-1$ grains (for any vertex $v$), i.e. the
darker it is, the closer to the stability threshold. Stable vertices have
greyscale colors, and unstable vertices have flashy colors. Table 1 presents
the color codes.

$deg(v)$ 0 grain 1 grain 2 grains 3 grains 4 grains 5 grains 6 grains 7 grains
3 4 6

Table 1:  Color codes of the pictures in this article, according to the degree
of the vertices (all our tilings have the same degree for all its
tiles/vertices).

##  1 Introduction

It all started in 1987 with the work of Bak, Tang and Wiesenfled [2].
Sandpiles have initially been introduced as number conserving cellular
automata on the two dimensional square grid, defined by the local toppling of
sand grains to neighboring cells, with statistics on chain reactions
presenting scale invariance typical of phase transition phenomena in physics,
the so called self-organized criticality [3, 34]. Soon after, Dhar realized
that sandpiles have a beautiful algebraic structure, which generalizes to
graphs and digraphs [10]. Since then, sandpiles have raised great interests
for their simple local definitions exhibiting complex global behaviors.
Following the work of Goles [23], numerous researches have been conducted on
one-dimensional models under sequential update mode [24, 25, 32, 50], parallel
update mode [13], and some variants such as symmetric [17, 20, 49, 51] or
Kadanoff rule [48].

On the two-dimensional side, the identity of the sandpile group on square
grids retains its mysteries, though the relaxation of hourglasses (toppling
from a single site) begin to reveal its structure through involved partial
differential equation developments [35, 36, 37, 44, 45].

The sandpile model on graphs is general enough to embed arbitrary computation,
as attested by its Turing-universality [27]. In order to analyse the apparent
complexity of the model in a formal framework, particular interest has been
raised on the computational complexity of predicting the dynamics of
sandpiles. Despite strong efforts started with Moore and Nilsson in [42] and
important contributions from Goles, the computational complexity of the
problem on two dimensional grids remains open between ${\mathsf{NC}}$ and
${\mathsf{P}}$. The one-dimensional model is in ${\mathsf{NC}}$ [40, 42], as
well as the Kadanoff rule [16, 18] and more general variants [19]. It is
${\mathsf{P}}$-complete from dimension three [42], and even allows for some
undecidable problems [6]. The (im)possibility of building crossing gates in
two dimensions is studied in [22, 43], and similar issues on closely related
threshold automata such as the majority rule have been fruitful [26, 28, 29,
30, 31, 41]. See [21] for a survey of the results on lattices.

Achieving isotropy in a cellular automaton is not a trivial matter, as cells
evolve on intrinsically anisotropic supports (grids). The authors of [9]
manage to build parabolas and circles in a two-dimensional cellular automaton
with $5^{13}$ states. Though optimizing the number of state was not an
objective of their work, it reflects the difficulty of the task. Approximative
(and simpler) approaches to build circles (isotropic diffusion) include:
probabilistic methods [39, 54, 55, 58], using a continuous state space [38],
or a large neighborhood to alleviate the anisotropy [14, 56, 59].

We give the definition of the sandpile model on general graphs in Section 2,
and present its algebraic structure at the heart of the experiments conducted
in this article. We also define the tilings considered as supports for the
sandpile dynamics: square, triangular and hexagonal grids, Penrose tilings
(kite-dart and rhombus) obtained by substitution, and Penrose tilings
(rhombus) obtained by the cut and project (multigrid) method. Identity
elements of the sandpile group on Penrose tilings are exposed in Section 3.
Our main contribution is in Section 4, where we study numerically the isotropy
observed during the stabilization process from the maximum stable
configuration plus the identity.

##  2 Model

Bak, Tang and Wiesenfeld first defined in [2] the sandpile model on two-
dimensional grids with von Neumann neighborhood. We present in Subsection 2.1
a general definition on any undirected multi graph, in Subsection 2.2 the
algebraic structure first revealed by Dhar in [10], and in Subsection 2.3 we
introduce various grids and more generally tilings on which the sandpile model
can be studied.

###  2.1 Sandpiles on graphs

Given a (connected) finite undirected multi graph $G=(V,E)$ with a
distinguished sink $s\in V$, let $\tilde{V}=V\setminus\\{s\\}$. A
configuration $c:\tilde{V}\to\mathbb{N}$ assigns a number of sand grains to
each non-sink vertex of $G$. The basic local evolution rule is the toppling at
some vertex $v\in\tilde{V}$, which may occur when $c(v)\geq deg(v)$, and
consists in vertex $v$ giving as many grains as it has edges to each of its
neighbors. When $c(v)\geq deg(v)$ we say that vertex $v$ is unstable, and a
configuration with no unstable vertex is called stable. Formally, let
$\Delta=D-A$ be the graph Laplacian of $G$, i.e. with $D$ the degree matrix
having $deg(v)$ on the diagonal and $A$ the adjacency matrix of $G$, and let
$\tilde{\Delta}$ be the reduced graph Laplacian obtained from $\Delta$ by
deleting the row and column corresponding to the sink $s$. With
$\tilde{\Delta}_{v}$ the row of $\tilde{\Delta}$ corresponding to vertex
$v\in\tilde{V}$, performing a toppling at $v$ corresponds to going from $c$ to
$c^{\prime}=c-\tilde{\Delta}_{v}$, which we denote
$c\overset{v}{\rightarrow}c^{\prime}$, or simply
$c\overset{}{\rightarrow}c^{\prime}$, and $\overset{}{\rightarrow}^{*}$ its
reflexive-transitive closure.

As such the system is non-deterministic on configuration space
$\mathbb{N}^{\tilde{V}}$, but it is straightforward to notice that any
configuration $c$ converges (because of the absorbing sink vertex $s$) to a
unique (because topplings commute) stable configuration, denoted
${c}^{\circ}$. This is the so called abelian property of sandpiles, or
fundamental theorem of sandpiles. The vector ${x}$ such that
$c-\tilde{\Delta}\cdot{x}={c}^{\circ}$ is called the shot vector or odometer
function of configuration $c$, it records how many times each vertex has
toppled during the stabilization of $c$. The fundamental theorem of sandpiles
furthermore states that the odometer function of any configuration is unique.

For the purpose of simulations that take some importance in the present work,
we define the deterministic synchronous dynamics as the evolution toppling
synchronously all unstable vertices, at each step. Formally, we define
$c\Rightarrow c^{\prime}$ with

$c^{\prime}=c-\sum_{\begin{subarray}{c}v\in\tilde{V}\\\ c(v)\geq
deg(v)\end{subarray}}\tilde{\Delta}_{v}.$

An example is given on Figure 1. For the sake of simplicity, in the following
graph will stand for finite undirected multi graph.

$\Rightarrow$ $\Rightarrow$ $\Rightarrow$ $\Rightarrow$ $\Rightarrow$

Figure 1:  Five steps of the sandpile model on a square grid of size $5\times
5$, corresponding to the graph with one vertex per square and north-east-
south-west adjacencies (all vertices have degree $4$). The sink $s$ is not
pictured, squares/vertices on the border have one edge connected to $s$, and
squares/vertices on the corners have two edges connected to $s$.

###  2.2 Abelian group structure

Given a graph $G=(V,E)$ with sink $s$, let
$\mathcal{C}=\mathbb{N}^{\tilde{V}}$ denote its set of configurations. The set
of stable configurations $\mathcal{C}_{\texttt{stab}}$ is naturally equipped
with the operation ${\,\oplus\,}$ defined as

$c{\,\oplus\,}c^{\prime}={(c+c^{\prime})}^{\circ}$

where $c+c^{\prime}$ is the componentwise addition of two configurations, i.e.
$(c+c^{\prime})(v)=c(v)+c^{\prime}(v)$ for all $v\in\tilde{V}$. It is
straightforward to notice that $(\mathcal{C}_{\texttt{stab}},{\,\oplus\,})$ is
a commutative monoid (closure, associativity, identity, commutativity), the
identity being the configuration $z$ such that $z(v)=0$ for all
$v\in\tilde{V}$.

Here comes the magics of sandpiles. The set of recurrent configurations

$\displaystyle\mathcal{C}_{\texttt{rec}}=$
$\displaystyle\\{c\in\mathcal{C}_{\texttt{stab}}\mid\forall
c^{\prime}\in\mathcal{C}:\exists
c^{\prime\prime}\in\mathcal{C}:c^{\prime}{\,\oplus\,}c^{\prime\prime}=c\\}$
$\displaystyle=$ $\displaystyle\\{c\in\mathcal{C}_{\texttt{stab}}\mid\forall
c^{\prime}\in\mathcal{C}_{\texttt{stab}}:\exists
c^{\prime\prime}\in\mathcal{C}_{\texttt{stab}}:c^{\prime}{\,\oplus\,}c^{\prime\prime}=c\\}$

corresponds to the intersection of ideals of $\mathcal{C}_{\texttt{stab}}$,
i.e.

$\mathcal{C}_{\texttt{rec}}=\bigcap_{\begin{subarray}{c}I\subseteq\mathcal{C}_{\texttt{stab}}\\\
I\text{ ideal of }\mathcal{C}_{\texttt{stab}}\end{subarray}}I$

with $I$ (right) ideal of $\mathcal{C}_{\texttt{stab}}$ if and only if
$I{\,\oplus\,}\mathcal{C}_{\texttt{stab}}\subseteq\mathcal{C}_{\texttt{stab}}$
where $I{\,\oplus\,}\mathcal{C}_{\texttt{stab}}=\\{i{\,\oplus\,}c\mid i\in
I\text{ and }c\in\mathcal{C}_{\texttt{stab}}\\}$. Moreover, it is a classical
result of algebra that the intersection of all ideals of a commutative
semigroup (closure, associativity, commutativity) gives an abelian group
(closure, associativity, identity, inverse, commutativity). So
$(\mathcal{C}_{\texttt{rec}},{\,\oplus\,})$ is an abelian group, called the
sandpile group on graph $G$ with sink $s$. Now remark that the configuration
$z$ containing no grain is (except on very restricted cases) not an element of
$\mathcal{C}_{\texttt{rec}}$, hence the111The identity element of the sandpile
group is unique. identity element $e\in\mathcal{C}_{\texttt{rec}}$ of the
sandpile group is a priori not obvious to construct, and it turns out that few
is know about its structure on numerous interesting cases, as presented on
Figure 2. It can be proven that

$e={(2m-{(2m)}^{\circ})}^{\circ}$ (1)

with $m$ the maximum stable configuration defined as $m(v)=deg(v)-1$ for all
$v\in\tilde{V}$, and $(2m)(v)=2\,c(v)$ for all $v\in\tilde{V}$. Indeed,
$2m-{(2m)}^{\circ}$ contains at least $deg(v)-1$ grains at each vertex
$v\in\tilde{V}$ hence its stabilization is recurrent, and furthermore it
corresponds to subtracting two configurations from the same equivalence class
according to relation $\leftrightarrow^{*}$ (the symmetric closure of
$\overset{}{\rightarrow}^{*}$), therefore to a configuration in the class of
the identity (see [11] for details). The equality follows since the identity
element of $(\mathcal{C}_{\texttt{rec}},{\,\oplus\,})$ is unique.

Figure 2:  Identity elements of the sandpile group on graphs corresponding to
the triangular grid of side length $50$ (left, $deg(v)=3$ for any vertex $v$),
the square grid of side length $50$ (center, $deg(v)=4$ for any vertex $v$),
the hexagonal grid of side length $14$ (right, $deg(v)=6$ for any vertex $v$).
All at the same scale.

###  2.3 Tilings

We will concentrate on finite tilings, but still introduce general
definitions. An infinite tiling $\mathcal{T}$ by $\tau$ is a covering of
$\mathbb{R}^{2}$ by finitely many polygonal tiles and their images by isometry
(translation, rotation, flip), i.e. copies of the tiles from $\tau$ cover the
plane without gaps nor overlaps. Let $\tau$ be a finite set of polygonal tiles
called a tile set, then $\mathcal{T}$ is a partition of $\mathbb{R}^{2}$ into
countably many isometries of the elements from $\tau$.

An infinite tiling $\mathcal{T}$ is periodic when it has a non-null
periodicity vector $\vec{p}\in\mathbb{R}^{2}$, such that
$\mathcal{T}+\vec{p}=\mathcal{T}$. A tile set $\tau$ is aperiodic when it
admits at least one infinite tiling (we say that $\tau$ tiles the plane), and
none are periodic. Aperiodicity is fundamentally related to the
uncomputability of the domino problem: given a (finite) tile set, does it tile
the plane? Let us quickly mention the seminal contributions of Wang [57],
Berger [5] and Robinson [52], along with the book Tilings and Patterns by
Grünbaum and Shephard [33].

In order to consider sandpiles on tilings, we explain now how to construct
finite undirected multi graphs with a distinguished sink. A finite tiling is
simply a subset of some infinite tiling $\mathcal{T}$. Remark that this
requires a finite tiling to be extensible into an infinite tiling, which will
be the case for all our finite tilings. Indeed, given some tile set $\tau$ we
will generate arbitrarily large finite tilings, which implies the existence of
an infinite tiling by compactness of the set of infinite tilings by $\tau$.
Two tiles are adjacent in $\mathcal{T}$ when they share an edge. All the
tilings we consider will see their adjacent tiles share full sides, i.e. no
tile will share a partial side or more than one side with another tile. These
are called edge-to-edge tilings. Given a finite tiling, each tile corresponds
to a vertex, plus an additional sink vertex corresponding to the outside face.
The tile to tile adjacencies are given by the edge-to-edge connections, and
the number of edges connecting a tile to the sink is equal to its number of
sides connected to the outside face, so that the degree of every tile is equal
to its number of sides. Tiles connected to the sink are said to be on the
border of the Tiling. Note that the sandpile dynamics is given by this
underlying graph, but tiles furthermore have coordinates in $\mathbb{R}^{2}$.

We now present the finite tilings considered in this article. The size of
tiles and the position of coordinate $(0,0)$ will be important for the
observations presented in Section 4. Grids are illustrated on Figure 2.

#### Square grids.

They correspond to the original two-dimensional sandpile model by Bak, Tang
and Wiesenfled [2]. Given size $n$, it basically consists in a $n\times n$
square grid with adjacencies given by von Neumann neighborhood (north, east,
south, west). Tiles on the borders have one edge connected to the sink, and
tiles on the corners have two edges connected to the sink. Each tile is a
square of side length $1$, and the finite tiling is centered with coordinate
$(0,0)$ in the middle of the grid: if $n$ is even then coordinate $(0,0)$
corresponds to four tiles corners; otherwise coordinate $(0,0)$ corresponds to
the center of a tile.

#### Triangular grids.

Given size $n$, it is made of equilateral triangles of side length $1$,
arranged up and down to form an equilateral triangle of side length $n$, where
the three outer sides of the tiling are made of $n$ triangular tiles. Tiles on
the border have one edge connected to the sink, and tiles on the corners have
two edges connected to the sink. The finite tiling is centered with coordinate
$(0,0)$ at the barycenter of its three corners.

#### Hexagonal grids.

Given size $n$, it is made of regular hexagons of side length $1$, arranged to
form an hexagonal grid (orientation is flat) with six sides each made of $n$
tiles. Tiles on the border have two edges connected to the sink, and tiles on
corners have three edges connected to the sink. The finite tiling is centered
with coordinate $(0,0)$ at the center of the central hexagon.

#### Penrose tilings.

Penrose developed in [46, 47] a series of elegant aperiodic tile sets. We
consider P2 (kite-dart), and P3 (rhombus) (tilings by P2 and P3 are mutually
locally derivable, see [1]). Penrose tilings may be obtained by
substitution222Note that in order to enforce aperiodicity in P2 and P3,
matching constraints should be added on tile edges, for example via notches,
but the finite tiling generation methods we employ do not require such
considerations. as described on Figure 4 for P2, and Figure 4 for P3. After
substituting, we rescale all tiles up by the substitution factor
$\phi=\frac{1+\sqrt{5}}{2}$, so that the tiles of the final tiling have the
exact same size as the base tiles (kite, dart, fat, thin). Figure 7 presents
three iterations of the substitution from a P2 Sun, Figure 7 presents three
iterations of the substitution from a P2 Star, and Figure 7 presents three
iterations of the substitution from a P3 Sun. Coordinate $(0,0)$ is at the
center of the initial Suns and Stars, and remains at the symmetry center of
tilings obtained by subsequent substitutions.

Figure 3:  Substitutions of P2 kite (left) and dart (right) tiles.

Figure 4:  Substitutions of P3 fat (left) and thin (right) rhombi tiles. Due
to the symmetry of P3 tiles, we highlight the origin point of each tile.

Figure 5:  Three iterations of the susbtitution from a P2 Sun, with identity
elements of the sandpile group pictured.

Figure 6:  Three iterations of the susbtitution from a P2 Star, with identity
elements of the sandpile group pictured.

Figure 7:  Three iterations of the susbtitution from a P3 Sun, with identity
elements of the sandpile group pictured.

Penrose P3 tilings may alternatively be generated by cut and project method:
consider the 5-dimensional plane $E$ spaned by vectors
$(\cos\frac{2k\pi}{5})_{0\leq k<5}$ and $(\sin\frac{2k\pi}{5})_{0\leq k<5}$,
passing through the point
$(\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5})$. The
orthogonal projection of the 5-dimensional grid lines (1-simplices) contained
in $E+[0,1]^{5}$ onto $E$ gives a tiling by P3 tiles. We implement it via the
dual multigrid method of de Bruijn [8], where one considers a 2-dimensional
space and five line families (our pentagrid) given by the intersections of $E$
with the five 5-dimensional hyperplanes $G_{i}=\\{x\in\mathbb{R}^{5}\mid
x\cdot e_{k}\in\mathbb{Z}\\}$ for $0\leq k<5$, where $e_{k}$ is the $k$-th
unit vector of $\mathbb{R}^{5}$. Pentagrid lines of a family are indexed by
the value of $x\cdot e_{k}$, and we denote $\ell^{k}_{i}$ the line from family
$0\leq k<5$ of index $i\in\mathbb{Z}$. Note that with this setting, no more
than two pentagrid lines intersect at a given position. The 2-dimensional
space is divided into polygonal cells delimited by pentagrid lines. Each cell
$p$ is labeled by a $5$-tuple of integers $(p_{k})_{0\leq k<5}$ such that cell
$p$ lies in between lines $\ell^{k}_{p_{k}}$ and $\ell^{k}_{p_{k}+1}$. To each
cell $p$ corresponds a point (a tile’s bound coordinate) in the 2-dimensional
space, at $\sum_{0\leq k<5}p_{k}(\cos\frac{2k\pi}{5},\sin\frac{2k\pi}{5})$. It
follows that to each intersecting pair of pentagrid lines corresponds a tile,
whose two edge orientations are given by the two line families. Details can be
found in [15]. We bound this process to a finite P3 cut and project tiling of
size $n$ by considering only tiles corresponding to points of intersection
lying in between $\ell^{i}_{-n}$ and $\ell^{i}_{n}$ for all families $0\leq
i<5$. The coordinate $(0,0)$ of this tiling is given by the cell labeled
$(0,0,0,0,0)$, and is therefore a tile bound (shared by multiple tiles).
Rhombus tiles have side length $1$ as on Figure 4. Figure 8 presents some P3
cut and project tilings.

Figure 8:  P3 cut and project tilings of sizes $1$, $2$, $3$ and $4$, with
identity elements of the sandpile group pictured.

##  3 Sandpile identity on Penrose tilings

We were curious to see what the identity element of the sandpile group would
look like on (finite) Penrose tilings, and it seems that no particular
structure appears333Well, this is a bit disappointing, but we think that it is
worth showing that it does not appear to be a fruitful research direction, or
maybe a more insightful reader would encounter something out there…. Examples
are presented on Figures 9, 10 and 11, JS-Sandpile computes them from Formula
1.

Figure 9:  Identity element of the sandpile group on the tiling obtained after
7 iterations of the substitution from a P2 Sun ($6710$ tiles).

Figure 10:  Identity element of the sandpile group on the tiling obtained
after 6 iterations of the substitution from a P3 Sun ($5415$ tiles).

Figure 11:  Identity element of the sandpile group on a P3 cut and project
tiling of size $10$ ($2440$ tiles).

We nevertheless discovered an interesting phenomenon on P3 cut and project
tilings, where the identity elements seem to display some stability. Indeed,
identity elements on successive sizes have a somewhat large central part of
the configuration in common. This obeservation is presented on Figure 12. We
may therefore conjecture a convergence of the sandpile identity on P3 cut and
project tilings: as the size $n$ increases, a larger part of the identity is
fixed (remains the same for all sizes $n^{\prime}\geq n$). Let us remark that
this phenomenon does not seem to take place on P2 nor P3 tilings obtained by
subsitution.

Figure 12:  Differences between successive identity elements of the sandpile
group on P3 cut and project tilings. From left to right, top to bottom, are
displayed the identity elements on P3 cut and project tilings of sizes $11$
($2900$ tiles), $12$ ($3520$ tiles), $13$ ($4155$ tiles), $14$ ($4790$ tiles),
$15$ ($5570$ tiles), $16$ ($6250$ tiles), $17$ ($7140$ tiles), $18$ ($8080$
tiles), $19$ ($8890$ tiles) and $20$ ($9940$ tiles), where the part of the
configuration which differs from the previous size (sizes $n$ is compared to
size $n-1$) are highlighted with redshifted colors.

##  4 Isotropic dynamics

On square, triangular, hexagonal grids and Penrose tilings, a very interesting
phenomenon appears444This also takes place on other tilings, outside the scope
of the present work. during the stabilization process

${(m+e)}^{\circ}=m$

where $m$ is the maximum stable configuration ($m(v)=deg(v)-1$ for all tile
$v$) and $e$ is the identity element of the sandpile group. Indeed, during the
last phase of the stabilization process leading back to $m$ (which is a
uniform configuration in these cases since the number of neighbors is
identical for all tiles), one can see the configuration $m$ reappear from the
outside (near the border) towards the center, outside a shrinking circular
shape, in a process step by step covering the whole configuration with tiles
containing $deg(v)-1$ grains. The first part of the stabilization process is
quite involved. Two illustrations are given on Figures 13 and 14.

Figure 13:  Stabilization of ${(m+e)}^{\circ}=m$ on a square grid of side
length $50$. From left to right, top to bottom, are displayed the
configurations every $100$ time steps (starting with $m+e$ at step $0$, ending
with time step $700$). The process converges to $m$ at step $707$.

Figure 14:  Stabilization of ${(m+e)}^{\circ}=m$ on $5$ iterations of the
substitution from a P2 Sun. From left to right, on top are displayed time
steps $0,50,100,150$ and at the bottom time steps $165,180,195,210$. The
process converges to $m$ at step $220$.

###  4.1 Roundness

In order to measure this phenomenon, we introduce the roundness as follows.
First, we partition a configuration into two parts: the outside part with all
tiles having $deg(v)-1$ grains connected to the border, and the inner part.
Given a tiling $G=(V,E)$ with sink $s$ and a configuration
$c:\tilde{V}\to\mathbb{N}$, let the maximum stable components,
$\texttt{MSC}(c)=\\{V_{1},V_{2},\dots,V_{k}\\}$, be the connected components
of tiles (from $\tilde{V}$) having $deg(v)-1$ grains. Then the outer tiles is
the set

$\texttt{outer}(c)=\bigcup\limits_{\begin{subarray}{c}V_{i}\in\texttt{MSC}(c)\\\
\exists v\in V_{i}:\\{v,s\\}\in E\end{subarray}}V_{i}$

and the inner tiles is the set
$\texttt{inner}(c)=\tilde{V}\setminus\texttt{outer}(c)$. From this partition
of the set of tiles, we are interested in the frontier between the outer and
inner tiles, and how close it is from a perfect circle. This has to do with
the coordinates of tiles in the Euclidean space $\mathbb{R}^{2}$, so let us
denote $\texttt{coord}(v)$ the set of coordinates of the bounds of some tile
$v\in\tilde{V}$, and $\texttt{body}(v)$ the subset of $\mathbb{R}^{2}$ covered
by the tile. Our convention regarding the sink $s$ is discussed below.
Regarding the circle, let us denote $\texttt{B}(r)$ the inside of the circle
of radius $r\in\mathbb{R}$ centered at coordinate $(0,0)$, i.e.

$\texttt{B}(r)=\\{(x,y)\in\mathbb{R}^{2}\mid\sqrt{x^{2}+y^{2}}\leq r\\}.$

We may therefore denote $\texttt{body}(v)\cap\texttt{B}(r)=\emptyset$ to state
that tile $v$ is entirely outside the circle of radius $r$ centered at
$(0,0)$, and $\texttt{body}(v)\subseteq\texttt{B}(r)$ (or equivalently
$\texttt{coord}(v)\subseteq\texttt{B}(r)$ since we deal with polygons and the
circle is convex) to state that tile $v$ is entirely inside the same circle.
We define the outer radius as the maximum scalar $r$ such that all outer tiles
are outside the circle of radius $r$ around the origin $(0,0)$ of the tiling,

$\underline{\texttt{r}}(c)=\max\\{r\in\mathbb{R}_{+}\mid\forall
v\in\texttt{outer}(c):\texttt{body}(v)\cap\texttt{B}(r)=\emptyset\\}$

and the inner radius as the minimum scalar $r$ such that all inner tiles are
inside the circle of radius $r$ around the origin $(0,0)$ of the tiling,

$\overline{\texttt{r}}(c)=\min\\{r\in\mathbb{R}_{+}\mid\forall
v\in\texttt{inner}(c):\texttt{body}(v)\subseteq\texttt{B}(r)\\}.$

In order to deal with the case $\texttt{outer}(c)=\emptyset$, which may for
example be the case on some configurations $(m+e)$, we add the convention that
the sink $s$ is an infinite tile covering all the space outside the tiling,
whose coordinates are the union of all bounds from tile edges adjacent to the
sink, and that it always belongs to $\texttt{outer}(c)$. As a consequence the
outer radius is upper bounded by the radius of the inscribed circle (inside
the finite tiling) with center at $(0,0)$. The case
$\texttt{inner}(c)=\emptyset$ is not problematic since the minimum defining
$\overline{\texttt{r}}(c)$ is taken on $\mathbb{R}_{+}$, and would therefore
equal $0$ as expected. The inner radius is upper bounded by the radius of the
circumscribed circle (outside the finite tiling) with center at $(0,0)$. The
roundness of a configuration $c$ is then measured as the difference between
the inner and outer radii,

$\texttt{r}(c)=\overline{\texttt{r}}(c)-\underline{\texttt{r}}(c).$

Note that $0\leq\texttt{r}(c)$, and we have $\texttt{r}(c)=0$ when the
frontier between inner and outer tiles is a perfect circle. Two examples of
roundness are given on Figure 15.

Figure 15:  Examples of roundness measures, circles of outer radius
$\underline{\texttt{r}}(c)$ in blue, inner radius $\overline{\texttt{r}}(c)$
in red. Left: square grid of side length $50$, after $500$ steps from $(m+e)$
(converges to $m$ in $707$ steps), $\texttt{r}(c)\approx 21.633-20.224=1.409$.
Right: $5$ iterations of the substition from a P2 Sun, after $150$ steps from
$(m+e)$ (converges to $m$ in $220$ steps), $\texttt{r}(c)\approx
16.662-14.729=1.933$.

###  4.2 Base roundnesses

Remark that since all our tiles are polygonal (with three to six sides), we
cannot expect to reach roundness $0$ (except when the stabilization process
has converged to $m$). To get some easy to interpret base values, we consider
the diameter of each tile, as the diameter of the circumscribed circle around
the tile (smallest radius of a circle having the tile entirely is its
interior). See Figure 16. The greatest diameter of some tiling’s tiles can be
interpreted as an upper bound on the best achievable roundness, for any
radius. Indeed, consider some tiling and a circle $C$ (of radius at most equal
to the inscribed radius), then all tiles having the center of their
circumscribed circle on or outside the circle $C$ can be set as outer tiles,
and all tiles having the center of their circumscribed circle inside the
circle $C$ can be set as inner tiles, which results in a roundness smaller
than twice the radius of the greatest tile’s radius, i.e. smaller than the
greatest tile’s diameter. We obtain the base roundnesses presented on Table 2.

Figure 16:  Diameters of all tiles, as the diameter of a circumscribed circle
(in purple). All tiles at the same scale.

Tiling Square grids Triangular grids Hexagonal grids P2 tilings P3 tilings
Base roundness $\underset{\approx 1.414}{\sqrt{2}}$ $\underset{\approx
1.155}{\frac{2\sqrt{3}}{3}}$ $2$ $2$ $\underset{\approx
1.902}{2\sin(\frac{2\pi}{5})}$

Table 2: Base roundnesses as the greatest diameter of some tiling’s tiles.

###  4.3 Plots

We now present plots of roundness measured during the stabilization process
${(m+e)}^{\circ}=m$ on the different tilings considered in this article. We
decompose the stabilization process from $m+e$ to $m$ into two phases:

  * $\circ$

phase 1: the dynamics is erratic,

  * $\circ$

phase 2: the set of inner tiles slowly shrinks, until reaching
$\texttt{inner}(m)=\emptyset$.

The beginning of phase 2 is defined as the first step such that the inner
radius is smaller or equal to the inscribed radius of the tiling (maximum
radius of a circle entirely inside the finite tiling, and centered at the
origin). Remark that at the beginning of phase 2, all tiles on the border of
the tiling are outer tiles, because the polygonal shape of any inner tile on
the border would otherwise lead to the inner radius being greater than the
inscribed radius of the tiling. An important observation is that, in all the
experiments presented in this article and performed during its preparation,
once in phase 2 with all border tiles as outer tiles, then all border tiles
remain outer tiles555They all remain stable with $deg(v)-1$ grains until
reaching $m$. Observe that any outer tile receiving some grain would topple,
and that toppling any outer tile would result in toppling the whole maximum
stable component it belongs to., and that the inner radius remains smaller or
equal to the inscribed radius. Two full examples of roundness plots are given
on Figures 18, 20, and their companion Figures 18, 20.

Figure 17:  Plot of the roundness during the stabilization process
${(m+e)}^{\circ}=m$, on the square grid of side length $50$. At each step we
plot the inner radius $\overline{\texttt{r}}(c)$ in red, outer radius
$\underline{\texttt{r}}(c)$ in blue, and roundness $\texttt{r}(c)$ in black.
Grid has one row per unit and one column per $10$ time steps. For example, at
step $3$, with $m+e\Rightarrow c^{1}\Rightarrow c^{2}\Rightarrow c^{3}$, we
observe that $\overline{\texttt{r}}(c^{3})=25\sqrt{2}\approx 35.355$ (the
radius of the circumscribed circle around the whole tiling) and that
$\underline{\texttt{r}}(c^{3})=25$ (the radius of the inscribed circle inside
the whole tiling), so that $\texttt{r}(c^{3})=25(\sqrt{2}-1)\approx 10.355$
(at step $0$ some outer tiles near the $x=0$ and $y=0$ axis give different
radii). Phase 2 begins at step $424$, and configuration $m$ is reached at step
$707$. The base roundness of this tiling (dashed) is $\sqrt{2}$. See also the
companion Figure 18.

Figure 18:  Configurations at steps $0$ (this is $m+e$), $423$ and $424$
during the stabilization process ${(m+e)}^{\circ}=m$, on the square grid of
side length $50$. Inner radii in red, outer radii in blue. One can observe the
phase transition occurring at step $424$: the inner radius becomes smaller or
equal to the inscribed radius of the tiling.

Figure 19:  Plot of the roundness during the stabilization process
${(m+e)}^{\circ}=m$, on the tiling obtained after $5$ iterations of the
substitution from a P2 Sun. At each step we plot the inner radius
$\overline{\texttt{r}}(c)$ in red, outer radius $\underline{\texttt{r}}(c)$ in
blue, and roundness $\texttt{r}(c)$ in black. Grid has one row per unit and
one column per $10$ time steps. For example, at step $0$ we observe that
$\overline{\texttt{r}}(m+e)\approx 22.940$ (the radius of the circumscribed
circle around the whole tiling) and that $\underline{\texttt{r}}(m+e)\approx
17.069$ (the radius of the inscribed circle inside the whole tiling), so that
$\texttt{r}(m+e)\approx 5.871$. Phase 2 begins at step $146$, and
configuration $m$ is reached at step $220$. The base roundness of this tiling
(dashed) is $2$. See also the companion Figure 20.

Figure 20:  Configurations at steps $0$ (this is $m+e$), $145$ and $146$
during the stabilization process ${(m+e)}^{\circ}=m$, on the tiling obtained
after $5$ iterations of the substitution from a P2 Sun. Inner radii in red,
outer radii in blue. One can observe the phase transition occurring at step
$146$: the inner radius becomes smaller or equal to the inscribed radius of
the tiling.

Experimental results, picturing only phase 2 of the stabilization process from
$m+e$ to $m$, are presented on Figure 21 and 22. Note that during all the
experiments we have performed in preparing this article, we have observed
similar behaviors for all other sizes.

Figure 21:  Plots of the roundness during phase 2 of the stabilization process
${(m+e)}^{\circ}=m$, on the square grid of side length $500$, triangular grid
of side length $500$, and hexagonal grid of side length $200$. Grids have one
row per unit and one column per $1000$ time steps. Two points are drawn every
$100$ time steps: the maximum and minimum roundness values $\texttt{r}(c)$
observed among these $100$ steps.

Figure 22:  Plots of the roundness during phase 2 of the stabilization process
${(m+e)}^{\circ}=m$, on the P2 Sun after $11$ iterations of the substitution,
P2 Star after $11$ iterations of the substitution, P3 Sun after $10$
iterations of the substitution, P3 cut and project tiling of size $100$. Grids
have one row per unit and one column per $1000$ time steps. Two points are
drawn every $100$ time steps: the maximum and minimum roundness values
$\texttt{r}(c)$ observed among these $100$ steps.

###  4.4 Computation times and data

The graphics were computed on our personal machines (standard laptops), with
simulation times and data presented on Table 3. Details on the implementation
(no parallelization) can be found at

https://github.com/huacayacauh/JS-Sandpile/wiki/Roundness.

Tiling Size/ iter. Tiles Identity (seconds) Roundness (seconds) Phase 2 begin
step Stabilization step Square grid $500$ $250000$ $825$ $1294$ $40413$
$70302$ Triangular grid $500$ $250000$ $1143$ $689$ $50163$ $69445$ Hexagonal
grid $200$ $120601$ $271$ $405$ $13252$ $24968$ P2 Sun subst. $11$ $327750$
$2931$ $1850$ $42724$ $67326$ P2 Star subst. $11$ $234410$ $944$ $640$ $27337$
$41918$ P3 Sun subst. $10$ $266860$ $2227$ $1570$ $42067$ $62333$ P3 cut &
project $100$ $249610$ $3145$ $2386$ $39210$ $66495$

Table 3:  Data regarding the computations generating the graphics presented on
Figures 21 and 22. Running times are given for the computation of identities
(from Formula 1) and roundness measures.

###  4.5 Analysis of roundness measures

Well, it appears clearly that the circular shapes observed during the
stabilization process ${(m+e)}^{\circ}=m$ are not asymptotically approaching
perfect circles. Indeed, in all experiments conducted on large tilings
(Figures 21 and 22), at the beginning of phase 2 the roundness is above the
base roundness, whereas the base roundness is an upper bound on the best
achievable roundness (see Subsection 4.2).

As the circular shapes shrink, it is normal to see the roundness decrease
until the value $0$ on configuration $m$ (all plots reach the minimum
roundness $0$ at stabilization step). On the other hand, it appears that
roundness and tiling size are correlated, meaning that as the size of the
tiling increases, configurations at the beginning of phase 2 are less round.
This is illustrated on Figure 23.

Figure 23:  Maximum roundness $\texttt{r}(c)$ measured during phase 2 of the
stabilization process ${(m+e)}^{\circ}=m$ on square grids of various side
lengths, illustrating the correlation between roundness and size.

Although they do not tend to perfect circles, we may admit that these shapes
are close to be, in regard of tiling sizes (Table 3). Increases of roundness
are visible to the naked eye on hexagonal and triangular grids which deviate
largely from a perfect circle, but are hardly noticeable without measurements
on other tilings. Figure 24 illustrates this with the frontier between outer
and inner tiles, on the configuration obtained at the beginning of phase 2
during the experiments from Subsection 4.3 (since these configurations are
quite large, we picture only the frontier in order to reduce the numerical
weight of the present document). If they are not perfect circles, then we may
naturally ask: what characterize these frontier shapes for each tiling?

Figure 24:  Frontier between inner and outer tiles on the configuration
obtained at the beginning of phase 2, as the set of edges shared by one inner
tile and one outer tile. Outer radius in blue, inner radius in red. Top:
square grid of side length $500$ at step $40413$ ($\texttt{r}(c)\approx
5.687$), triangular grid of side length $500$ at step $50163$
($\texttt{r}(c)\approx 7.755$), hexagonal grid of side length $200$ at step
$13252$ ($\texttt{r}(c)\approx 14.470$). Middle: P2 Sun after $11$ iterations
of the subsitution at step $42724$ ($\texttt{r}(c)\approx 4.246$), P2 Star
after $11$ iterations of the subsitution at step $27337$
($\texttt{r}(c)\approx 4.064$). Bottom: P3 Sun after $10$ iterations of the
subsitution at step $42067$ ($\texttt{r}(c)\approx 3.220$), P3 Sun cut and
project tiling of size $100$ at step $39210$ ($\texttt{r}(c)\approx 2.452$).

Finally, we see on Figure 24 (first line) that the frontier on grids
(respectively square, triangular and hexagonal) reflect the shape on the
border of the tilings, somehow “rounded”. Indeed, it appears that the number
of corners of the grids are equal to the number of parts where the frontier
deviates from a circle (the number of times it goes from the inner radius to
the outer radius). Theses symmetries are expected, they come from the symmetry
of the grids and therefore the symmetries of the dynamics. A natural attempt
would be to experiment the roundness of a square grid cropped to a circle, in
order to remove the effect of the border’s shape. Despite the fact that the
difference between the circumscribed and inscribed radius is smaller than the
base roundness, this does not lead to significantly smaller roundness
measurements (frontiers are not closer to perfect circles) during phase 2 of
the stabilization process ${(m+e)}^{\circ}=m$, as shown on Figure 25. It feels
that the frontier reflects the anisotropy of the square grid itself rather
than that of its border’s shape.

Figure 25:  Left: maximum roundness $\texttt{r}(c)$ measured during phase 2 of
the stabilization process ${(m+e)}^{\circ}=m$ on square grids cropped to
circles of various radii. Right: frontier between inner and outer tiles on the
configuration obtained at the beginning of phase 2 (step $32131$,
$\texttt{r}(c)\approx 4.073$) on a square grid cropped to a circle of radius
$250$ (the border of the tiling is pictured in green).

##  5 Conclusions and perspectives

The experiments presented in this article are reproducible with JS-Sandpile
(links in Preamble). The software implements no parallelization mechanism,
which would allow to perform larger simulations (e.g. using GPU). Nevertheless
we believe that this would not lead to qualitatively different observations.

We have presented some identity elements on Penrose tilings, revealing no
obvious structure related to these famous aperiodic tilings. Identities are
highly sensitive to the shape of the tiling, and it may be the case that other
finite croppings of (infinite) Penrose tilings lead to different observations.
We tried to build the most “natural” finite Penrose tilings: Suns and Stars
obtained by substitution, along with cut and project from a $5$-dimensional
hypercube.

The apparent isotropy, observed during the stabilization of the maximum stable
configuration plus the identity element of the sandpile group, has been
measured through the notion of roundness. Experiments revealed that these
frontiers are actually not approaching perfect circles on any tiling under
consideration. Two further directions may be investigated. First, if these
shapes are not perfect circles, then how to characterize them on each tiling
(especially on grids)? Second, which tilings would lead to frontiers
approaching perfect circles (if any)? The modest attempt to crop a square grid
to a circle fails, suggesting that this may not be easy to achieve on tilings
which are intrinsically anisotropic at the tile level.

Penrose tilings exhibit frontiers closer to perfect circles, though they also
deviate significantly from their best achievable roundness. It is not very
surprising to observe on Penrose tilings a behavior similar to regular
lattices, as they are quasi-periodic, and in some sense the “most regular”
aperiodic tilings. Let us open the large perspective of considering other
tilings, for example tilings with higher quasi-periodicity functions (above
affine [4, 7, 12]666The difficulty may be to find constructions from non Wang
tiles, because Wang tiles would lead to square grids for the sandpile model to
play on (as we remove tile decorations).), Cayley graphs [53], or hyperbolic
planes (e.g. with Poincaré disk models).

Finally, could the stabilization process ${(m+e)}^{\circ}=m$ shed some light
on the enigmatic identity elements on grids? This is really out of reach for
our present knowledge, but we hope that the recent progresses of Levine et al.
[35, 36, 37, 44, 45] are breaking some scientific locks in the domain of
sandpiles.

## Acknowledgments

The authors are thankful to Valentin Darrigo for his contributions to JS-
Sandpile, to Victor Poupet for stressing that the apparent isotropy of the
${(m+e)}^{\circ}$ process is surprising at the occasion of a talk given by KP
during AUTOMATA’2014 in Himeji, to Thomas Fernique for sharing his expertise
(and code!) regarding the cut and project method, and to Christophe Papazian
for useful comments on quasi-periodicity functions.

The work of JF was conducted while a Master student at Aix-Marseille
Université, doing an internship at the LIS laboratory (UMR 7020), both in
Marseille, France. The work of KP was funded mainly by his salary as a French
State agent and therefore by French taxpayers’ taxes, affiliated to Aix-
Marseille Univ., Univ. de Toulon, CNRS, LIS, UMR 7020, Marseille, France and
Univ. Côte d’Azur, CNRS, I3S, UMR 7271, Sophia Antipolis, France. Secondary
financial support came from ANR-18-CE40-0002 FANs project, ECOS-Sud C16E01
project, and STIC AmSud CoDANet 19-STIC-03 (Campus France 43478PD) project.

## References

  * [1] M. Baake, M. Scholottmann, and P. D. Jarvis.  Quasiperiodic tilings with tenfold symmetry and equivalence with respect to local derivability.  Journal of Physics A: Mathematical and General, 24(19):4637–4654, 1991.
  * [2] P. Bak, C. Tang, and K. Wiesenfeld.  Self-organized criticality: An explanation of the 1/f noise.  Physical Review Letters, 59:381–384, 1987.
  * [3] P. Bak, C. Tang, and K. Wiesenfeld.  Self-organized criticality.  Physical Review A, 38(1):364–374, 1988.
  * [4] A. Ballier and E. Jeandel.  Computing (or not) Quasi-periodicity Functions of Tilings.  In Journées Automates Cellulaires 2010, pages 54–64, 2010\.
  * [5] R. Berger.  The undecidability of the domino problem.  Memoirs of the American Mathematical Society, 66, 1966.
  * [6] H. Cairns.  Some halting problems for abelian sandpiles are undecidable in dimension three.  SIAM Journal on Discrete Mathematics, 32(4):2636–2666, 2018.
  * [7] J. Cervelle and B. Durand.  Tilings: recursivity and regularity.  Theoretical Computer Science, 310(1):469–477, 2004.
  * [8] N. G. de Bruijn.  Algebraic theory of penrose’s non-periodic tilings of the plane. ii.  Indagationes Mathematicae, 84(1):53–66, 1981.
  * [9] M. Delorme, J. Mazoyer, and L. Tougne.  Discrete parabolas and circles on 2D cellular automata.  Theoretical Computer Science, 218(2):347–417, 1999.
  * [10] D. Dhar.  Self-organized critical state of sandpile automaton models.  Physical Review Letters, 64:1613–1616, 1990.
  * [11] D. Dhar, P. Ruelle, S. Sen, and D. N. Verma.  Algebraic aspects of abelian sandpile models.  Journal of Physics A, 28(4):805–831, 1995.
  * [12] B. Durand.  Tilings and quasiperiodicity.  Theoretical Computer Science, 221(1):61–75, 1999.
  * [13] J. O. Durand-Lose.  Parallel transient time of one-dimensional sand pile.  Theoretical Computer Science, 205(1-2):183–193, 1998.
  * [14] V. G. Fast and I. R. Efimov.  Stability of vortex rotation in an excitable cellular medium.  Physica D: Nonlinear Phenomena, 49(1):75–81, 1991.
  * [15] T. Fernique.  Pavages, Fractions Continues et Géométrie Discrète.  PhD thesis, Université Montpellier II, 2007.
  * [16] E. Formenti, E. Goles, and B. Martin.  Computational complexity of avalanches in the Kadanoff sandpile model.  Fundamenta Informaticae, 115(1):107–124, 2012.
  * [17] E. Formenti, B. Masson, and T. Pisokas.  Advances in symmetric sandpiles.  Fundamenta Informaticae, 76(1-2):91–112, 2007.
  * [18] E. Formenti, K. Perrot, and E. Rémila.  Computational complexity of the avalanche problem on one dimensional Kadanoff sandpiles.  In Proceedings of AUTOMATA’2014, volume 8996 of LNCS, pages 21–30, 2014.
  * [19] E. Formenti, K. Perrot, and E. Rémila.  Computational complexity of the avalanche problem for one dimensional decreasing sandpiles.  Journal of Cellular Automata, 13:215–228, 2018.
  * [20] E. Formenti, V. T. Pham, H. D. Phan, and T. H. Tran.  Fixed-point forms of the parallel symmetric sandpile model.  Theoretical Computer Science, 533:1–14, 2014.
  * [21] Enrico Formenti and Kévin Perrot.  How Hard is it to Predict Sandpiles on Lattices? A Survey.  Fondamenta Informaticae, 171:189–219, 2019.
  * [22] A. Gajardo and E. Goles.  Crossing information in two-dimensional sandpiles.  Theoretical Computer Science, 369(1-3):463–469, 2006.
  * [23] E. Goles.  Sand pile automata.  Annales de l’institut Henri Poincaré (A) Physique théorique, 56(1):75–90, 1992.
  * [24] E. Goles and M. Kiwi.  Games on line graphs and sand piles.  Theoretical Computer Science, 115(2):321–349, 1993.
  * [25] E. Goles, M. Latapy, C. Magnien, M. Morvan, and H. D. Phan.  Sandpile models and lattices: a comprehensive survey.  Theoretical Computer Science, 322(2):383–407, 2004.
  * [26] E. Goles, D. Maldonado, P. Montealegre, and N. Ollinger.  On the computational complexity of the freezing non-strict majority automata.  In Proceedings of AUTOMATA’2017, pages 109–119, 2017.
  * [27] E. Goles and M. Margenstern.  Universality of the chip-firing game.  Theoretical Computer Science, 172(1-2):121–134, 1997.
  * [28] E. Goles and P. Montealegre.  Computational complexity of threshold automata networks under different updating schemes.  Theoretical Computer Science, 559:3–19, 2014.
  * [29] E. Goles and P. Montealegre.  A fast parallel algorithm for the robust prediction of the two-dimensional strict majority automaton.  In Proceedings of ACRI’2016, pages 166–175, 2016.
  * [30] E. Goles, P. Montealegre, K. Perrot, and G. Theyssier.  On the complexity of two-dimensional signed majority cellular automata.  Journal of Computer and System Sciences, 91:1–32, 2017.
  * [31] E. Goles, P. Montealegre-Barba, and I. Todinca.  The complexity of the bootstraping percolation and other problems.  Theoretical Computer Science, 504:73–82, 2013.
  * [32] E. Goles, M. Morvan, and H. D. Phan.  Sandpiles and order structure of integer partitions.  Discrete Applied Mathematics, 117(1–3):51–64, 2002.
  * [33] B. Grünbaum and G. C. Shephard.  Tilings and Patterns.  W. H. Freeman & Co., 1986.
  * [34] L. P. Kadanoff, S. R. Nagel, L. Wu, and S. Zhou.  Scaling and universality in avalanches.  Physical Review A, 39(12):6524–6537, 1989.
  * [35] L. Levine, W. Pegden, and C. K. Smart.  Apollonian structure in the abelian sandpile.  Geometric and Functional Analysis, 26:306–336, 2016.
  * [36] L. Levine, W. Pegden, and C. K. Smart.  The apollonian structure of integer superharmonic matrices.  Annals of Mathematics, 186(1):1–67, 2017.
  * [37] L. Levine and Y. Peres.  Laplacian growth, sandpiles, and scaling limits.  Bulletin of the American Mathematical Society, 54(3):355–382, 2017\.
  * [38] M. Marek.  Grid anisotropy reduction for simulation of growth processes with cellular automaton.  Physica D: Nonlinear Phenomena, 253:73–84, 2013.
  * [39] M. Markus and B. Hess.  Isotropic cellular automaton for modelling excitable media.  Nature, 347:56–58, 1990.
  * [40] P. B. Miltersen.  The computational complexity of one-dimensional sandpiles.  In Proceedings of CiE’2005, pages 342–348, 2005.
  * [41] C. Moore.  Majority-vote cellular automata, Ising dynamics, and P-completeness.
  * [42] C. Moore and M. Nilsson.  The computational complexity of sandpiles.  Journal of Statistical Physics, 96:205–224, 1999.  10.1023/A:1004524500416.
  * [43] V.-H. Nguyen and K. Perrot.  Any shape can ultimately cross information on two-dimensional abelian sandpile models.  In Proceedings of AUTOMATA’2018, volume 10875 of LNCS, pages 127–142, 2018.
  * [44] W. Pegden and C. K. Smart.  Convergence of the abelian sandpile.  Duke Mathematical Journal, 162(4):627–642, 03 2013.
  * [45] W. Pegden and C. K. Smart.  Stability of patterns in the abelian sandpile.  Annales Henri Poincaré, 21:1383–1399, 2020.
  * [46] R. Penrose.  The role of aesthetics in pure and applied mathematical research.  Bulletin of the Institute of Mathematics and Its Applications, 10(2):266–271, 1974.
  * [47] R. Penrose.  Pentaplexity: A class of non-periodic tilings of the plane.  The Mathematical Intelligencer 2, pages 32––37, 1979.
  * [48] K. Perrot.  Les piles de sable Kadanoff.  PhD thesis, École normale supérieure de Lyon, 2013.
  * [49] K. Perrot, H. D. Phan, and V. T. Pham.  On the set of Fixed Points of the Parallel Symmetric Sand Pile Model.  In Proceedings AUTOMATA’2011, DMTCS, pages 17–28. Open Publishing Association, 2011.
  * [50] H. D. Phan.  Structures ordonnées et dynamiques de piles de sable.  PhD thesis, Université Paris 7, 1999.
  * [51] H. D. Phan.  Two sided sand piles model and unimodal sequences.  ITA, 42(3):631–646, 2008.
  * [52] R. M. Robinson.  Undecidability and nonperiodicity for tilings of the plane.  Inventiones mathematicae, 12:177–209, 1971.
  * [53] Zs. Roka.  Automates cellulaires sur graphes de Cayley.  PhD thesis, École Normale Supérieure de Lyon, 1994.
  * [54] H. E. Schepers and M. Markus.  Two types of performance of an isotropic cellular automaton: stationary (Turing) patterns and spiral waves.  Physica A: Statistical Mechanics and its Applications, 188(1):337–343, 1992.
  * [55] B. Schönfisch.  Anisotropy in cellular automata.  Biosystems, 41(1):29–41, 1997.
  * [56] G. Ch. Sirakoulis, I. Karafyllidis, and A. Thanailakis.  A cellular automaton for the propagation of circular fronts and its applications.  Engineering Applications of Artificial Intelligence, 18(6):731–744, 2005.
  * [57] H. Wang.  Proving theorems by pattern recognition — II.  The Bell System Technical Journal, 40:1–41, 1961.
  * [58] J. R. Weimar.  Cellular automata for reaction-diffusion systems.  Parallel Computing, 23(11):1699–1715, 1997.
  * [59] J. R. Weimar, J. J. Tyson, and L. T. Watson.  Diffusion and wave propagation in cellular automaton models of excitable media.  Physica D: Nonlinear Phenomena, 55(3):309–327, 1992.