Sandpile toppling on Penrose tilings: identity and isotropic dynamics
SSandpile toppling on Penrose tilings:identity and isotropic dynamics
J´er´emy Fersula , Camille Noˆus , and K´evin Perrot Universit´e publique Cogitamus laboratory
Abstract
We present experiments of sandpiles on grids (square, triangular, hexagonal) andPenrose tilings. The challenging part is to program such simulator; and our javacriptcode is available online, ready to play! We first present some identity elements ofthe sandpile group on these aperiodic structures, and then study the stabilization ofthe maximum stable configuration plus the identity, which lets a surprising circularshape appear. Roundness measurements reveal that the shapes are not approachingperfect circles, though they are close to be. We compare numerically this almostisotropic dynamical phenomenon on various tilings.
Preamble
The experiments presented in this paper were conducted using
JS-Sandpile , a javascriptsandpile 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 ) − v ), i.e. the darker it is, the closer tothe stability threshold. Stable vertices have greyscale colors, and unstable vertices haveflashy colors. Table 1 presents the color codes. deg ( v ) 0 grain grain grains grains grains grains grains grains a r X i v : . [ n li n . C G ] J un Introduction
It all started in 1987 with the work of Bak, Tang and Wiesenfled [2]. Sandpiles haveinitially been introduced as number conserving cellular automata on the two dimen-sional square grid, defined by the local toppling of sand grains to neighboring cells,with statistics on chain reactions presenting scale invariance typical of phase transitionphenomena in physics, the so called self-organized criticality [3, 34]. Soon after, Dharrealized that sandpiles have a beautiful algebraic structure, which generalizes to graphsand digraphs [10]. Since then, sandpiles have raised great interests for their simple lo-cal definitions exhibiting complex global behaviors. Following the work of Goles [23],numerous researches have been conducted on one-dimensional models under sequentialupdate mode [24, 25, 32, 50], parallel update mode [13], and some variants such assymmetric [17, 20, 49, 51] or Kadanoff rule [48].On the two-dimensional side, the identity of the sandpile group on square grids retainsits mysteries, though the relaxation of hourglasses (toppling from a single site) begin toreveal 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, asattested by its Turing-universality [27]. In order to analyse the apparent complexity ofthe model in a formal framework, particular interest has been raised on the computa-tional complexity of predicting the dynamics of sandpiles. Despite strong efforts startedwith Moore and Nilsson in [42] and important contributions from Goles, the computa-tional complexity of the problem on two dimensional grids remains open between NC and P . The one-dimensional model is in NC [40, 42], as well as the Kadanoff rule [16, 18]and more general variants [19]. It is P -complete from dimension three [42], and evenallows for some undecidable problems [6]. The (im)possibility of building crossing gatesin two dimensions is studied in [22, 43], and similar issues on closely related thresholdautomata 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 onintrinsically anisotropic supports (grids). The authors of [9] manage to build parabolasand circles in a two-dimensional cellular automaton with 5 states. Though optimizingthe 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 largeneighborhood to alleviate the anisotropy [14, 56, 59].We give the definition of the sandpile model on general graphs in Section 2, andpresent 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 bysubstitution, and Penrose tilings (rhombus) obtained by the cut and project (multigrid)method. Identity elements of the sandpile group on Penrose tilings are exposed inSection 3. Our main contribution is in Section 4, where we study numerically the isotropy2bserved during the stabilization process from the maximum stable configuration plusthe identity. Bak, Tang and Wiesenfeld first defined in [2] the sandpile model on two-dimensional gridswith von Neumann neighborhood. We present in Subsection 2.1 a general definition onany undirected multi graph, in Subsection 2.2 the algebraic structure first revealed byDhar in [10], and in Subsection 2.3 we introduce various grids and more generally tilingson which the sandpile model can be studied.
Given a (connected) finite undirected multi graph G = ( V, E ) with a distinguished sink s ∈ V , let ˜ V = V \ { s } . A configuration c : ˜ V → N assigns a number of sand grains toeach non-sink vertex of G . The basic local evolution rule is the toppling at some vertex v ∈ ˜ V , which may occur when c ( v ) ≥ deg ( v ), and consists in vertex v giving as manygrains as it has edges to each of its neighbors. When c ( v ) ≥ deg ( v ) we say that vertex v is unstable , and a configuration with no unstable vertex is called stable . Formally, let∆ = D − A be the graph Laplacian of G , i.e. with D the degree matrix having deg ( v ) onthe diagonal and A the adjacency matrix of G , and let ˜∆ be the reduced graph Laplacian obtained from ∆ by deleting the row and column corresponding to the sink s . With ˜∆ v the row of ˜∆ corresponding to vertex v ∈ ˜ V , performing a toppling at v correspondsto going from c to c (cid:48) = c − ˜∆ v , which we denote c v → c (cid:48) , or simply c → c (cid:48) , and → ∗ itsreflexive-transitive closure.As such the system is non-deterministic on configuration space N ˜ V , but it is straight-forward to notice that any configuration c converges (because of the absorbing sink vertex s ) to a unique (because topplings commute) stable configuration, denoted c ◦ . This is theso called abelian property of sandpiles, or fundamental theorem of sandpiles. The vector x such that c − ˜∆ · x = c ◦ 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 . Thefundamental theorem of sandpiles furthermore states that the odometer function of anyconfiguration is unique.For the purpose of simulations that take some importance in the present work, wedefine the deterministic synchronous dynamics as the evolution toppling synchronouslyall unstable vertices, at each step. Formally, we define c ⇒ c (cid:48) with c (cid:48) = c − (cid:88) v ∈ ˜ Vc ( v ) ≥ deg ( v ) ˜∆ v . An example is given on Figure 1. For the sake of simplicity, in the following graph willstand for finite undirected multi graph . 3 ⇒ ⇒ ⇒ ⇒
Figure 1: Five steps of the sandpile model on a square grid of size 5 ×
5, corresponding tothe graph with one vertex per square and north-east-south-west adjacencies (all verticeshave degree 4). The sink s is not pictured, squares/vertices on the border have one edgeconnected to s , and squares/vertices on the corners have two edges connected to s . Given a graph G = ( V, E ) with sink s , let C = N ˜ V denote its set of configurations. Theset of stable configurations C stab is naturally equipped with the operation ⊕ defined as c ⊕ c (cid:48) = ( c + c (cid:48) ) ◦ where c + c (cid:48) is the componentwise addition of two configurations, i.e. ( c + c (cid:48) )( v ) = c ( v ) + c (cid:48) ( v ) for all v ∈ ˜ V . It is straightforward to notice that ( C stab , ⊕ ) is a commu-tative monoid (closure, associativity, identity, commutativity), the identity being theconfiguration z such that z ( v ) = 0 for all v ∈ ˜ V .Here comes the magics of sandpiles. The set of recurrent configurations C rec = { c ∈ C stab | ∀ c (cid:48) ∈ C : ∃ c (cid:48)(cid:48) ∈ C : c (cid:48) ⊕ c (cid:48)(cid:48) = c } = { c ∈ C stab | ∀ c (cid:48) ∈ C stab : ∃ c (cid:48)(cid:48) ∈ C stab : c (cid:48) ⊕ c (cid:48)(cid:48) = c } corresponds to the intersection of ideals of C stab , i.e. C rec = (cid:92) I ⊆C stab I ideal of C stab I with I (right) ideal of C stab if and only if I ⊕ C stab ⊆ C stab where I ⊕ C stab = { i ⊕ c | i ∈ I and c ∈ C stab } . Moreover, it is a classical result of algebra that the intersectionof all ideals of a commutative semigroup (closure, associativity, commutativity) gives an abelian group (closure, associativity, identity, inverse, commutativity). So ( C rec , ⊕ ) isan abelian group, called the sandpile group on graph G with sink s . Now remark that theconfiguration z containing no grain is (except on very restricted cases) not an element of C rec , hence the identity element e ∈ C rec of the sandpile group is a priori not obvious toconstruct, and it turns out that few is know about its structure on numerous interestingcases, as presented on Figure 2. It can be proven that e = (2 m − (2 m ) ◦ ) ◦ (1)with m the maximum stable configuration defined as m ( v ) = deg ( v ) − v ∈ ˜ V , and(2 m )( v ) = 2 c ( v ) for all v ∈ ˜ V . Indeed, 2 m − (2 m ) ◦ contains at least deg ( v ) − The identity element of the sandpile group is unique. deg ( v ) = 3 for any vertex v ), the square grid ofside 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.each vertex v ∈ ˜ V hence its stabilization is recurrent, and furthermore it corresponds tosubtracting two configurations from the same equivalence class according to relation ↔ ∗ (the symmetric closure of → ∗ ), therefore to a configuration in the class of the identity(see [11] for details). The equality follows since the identity element of ( C rec , ⊕ ) isunique. We will concentrate on finite tilings, but still introduce general definitions. An infinitetiling T by τ is a covering of R by finitely many polygonal tiles and their images byisometry (translation, rotation, flip), i.e. copies of the tiles from τ cover the planewithout gaps nor overlaps. Let τ be a finite set of polygonal tiles called a tile set , then T is a partition of R into countably many isometries of the elements from τ .An infinite tiling T is periodic when it has a non-null periodicity vector (cid:126)p ∈ R ,such that T + (cid:126)p = T . A tile set τ is aperiodic when it admits at least one infinitetiling (we say that τ tiles the plane ), and none are periodic. Aperiodicity is fundamen-tally 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¨unbaumand Shephard [33].In order to consider sandpiles on tilings, we explain now how to construct finiteundirected multi graphs with a distinguished sink. A finite tiling is simply a subset ofsome infinite tiling T . Remark that this requires a finite tiling to be extensible into aninfinite tiling, which will be the case for all our finite tilings. Indeed, given some tileset τ we will generate arbitrarily large finite tilings, which implies the existence of aninfinite tiling by compactness of the set of infinite tilings by τ . Two tiles are adjacent in T when they share an edge. All the tilings we consider will see their adjacent tilesshare full sides, i.e. no tile will share a partial side or more than one side with another5ile. These are called edge-to-edge tilings. Given a finite tiling, each tile correspondsto a vertex, plus an additional sink vertex corresponding to the outside face . The tileto tile adjacencies are given by the edge-to-edge connections, and the number of edgesconnecting 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 thesink are said to be on the border of the Tiling. Note that the sandpile dynamics is givenby this underlying graph, but tiles furthermore have coordinates in R .We now present the finite tilings considered in this article. The size of tiles and theposition 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 byBak, Tang and Wiesenfled [2]. Given size n , it basically consists in a n × n square gridwith adjacencies given by von Neumann neighborhood (north, east, south, west). Tileson the borders have one edge connected to the sink, and tiles on the corners have twoedges connected to the sink. Each tile is a square of side length 1, and the finite tilingis 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 thecenter 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 threeouter sides of the tiling are made of n triangular tiles. Tiles on the border have one edgeconnected 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, ar-ranged 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 threeedges connected to the sink. The finite tiling is centered with coordinate (0 ,
0) at thecenter 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 locallyderivable , see [1]). Penrose tilings may be obtained by substitution as described onFigure 3 for P2, and Figure 4 for P3. After substituting, we rescale all tiles up bythe substitution factor φ = √ , so that the tiles of the final tiling have the exactsame size as the base tiles (kite, dart, fat, thin). Figure 5 presents three iterations ofthe substitution from a P2 Sun , Figure 6 presents three iterations of the substitutionfrom a
P2 Star , and Figure 7 presents three iterations of the substitution from a P3 Note that in order to enforce aperiodicity in P2 and P3, matching constraints should be added ontile edges, for example via notches, but the finite tiling generation methods we employ do not requiresuch considerations. φ π π φ φ π π 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 thesymmetry of P3 tiles, we highlight the origin point of each tile.
Sun . Coordinate (0 ,
0) is at the center of the initial Suns and Stars, and remains at thesymmetry center of tilings obtained by subsequent substitutions.Penrose P3 tilings may alternatively be generated by cut and project method: con-sider the 5-dimensional plane E spaned by vectors (cos kπ ) ≤ k< and (sin kπ ) ≤ k< ,passing through the point ( , , , , ). The orthogonal projection of the 5-dimensionalgrid lines (1-simplices) contained in E + [0 , onto E gives a tiling by P3 tiles. Weimplement 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 ∈ R | x · e k ∈ Z } for 0 ≤ k <
5, where e k is the k -th unit vector of R . Pentagrid lines of a family are indexed by the valueof x · e k , and we denote (cid:96) ki the line from family 0 ≤ k < i ∈ Z . Note thatwith this setting, no more than two pentagrid lines intersect at a given position. The2-dimensional space is divided into polygonal cells delimited by pentagrid lines. Eachcell p is labeled by a 5-tuple of integers ( p k ) ≤ k< such that cell p lies in between lines (cid:96) kp k and (cid:96) kp k +1 . To each cell p corresponds a point (a tile’s bound coordinate) in the2-dimensional space, at (cid:80) ≤ k< p k (cos kπ , sin kπ ). It follows that to each intersectingpair of pentagrid lines corresponds a tile, whose two edge orientations are given by thetwo line families. Details can be found in [15]. We bound this process to a finite P3 cutand project tiling of size n by considering only tiles corresponding to points of intersec-tion lying in between (cid:96) i − n and (cid:96) in for all families 0 ≤ i <
5. The coordinate (0 ,
0) of thistiling is given by the cell labeled (0 , , , , P2 Sun , with identity elements ofthe sandpile group pictured.Figure 6: Three iterations of the susbtitution from a
P2 Star , with identity elementsof the sandpile group pictured.Figure 7: Three iterations of the susbtitution from a
P3 Sun , with identity elements ofthe sandpile group pictured. 8igure 8: P3 cut and project tilings of sizes 1, 2, 3 and 4, with identity elements of thesandpile group pictured.
We were curious to see what the identity element of the sandpile group would look likeon (finite) Penrose tilings, and it seems that no particular structure appears . Examplesare presented on Figures 9, 10 and 11, JS-Sandpile computes them from Formula 1.We nevertheless discovered an interesting phenomenon on P3 cut and project tilings,where the identity elements seem to display some stability. Indeed, identity elements onsuccessive sizes have a somewhat large central part of the configuration in common. Thisobeservation is presented on Figure 12. We may therefore conjecture a convergence ofthe sandpile identity on P3 cut and project tilings: as the size n increases, a larger partof the identity is fixed (remains the same for all sizes n (cid:48) ≥ n ). Let us remark that thisphenomenon does not seem to take place on P2 nor P3 tilings obtained by subsitution. On square, triangular, hexagonal grids and Penrose tilings, a very interesting phe-nomenon appears during the stabilization process( m + e ) ◦ = m where m is the maximum stable configuration ( m ( v ) = deg ( v ) − v ) and e is the identity element of the sandpile group. Indeed, during the last phase of thestabilization process leading back to m (which is a uniform configuration in these casessince 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 shrinkingcircular shape, in a process step by step covering the whole configuration with tiles Well, this is a bit disappointing, but we think that it is worth showing that it does not appear to bea fruitful research direction, or maybe a more insightful reader would encounter something out there. . . This also takes place on other tilings, outside the scope of the present work. n is compared to size n −
1) are highlighted with redshifted colors.13igure 13: Stabilization of ( m + e ) ◦ = m on a square grid of side length 50. From leftto right, top to bottom, are displayed the configurations every 100 time steps (startingwith m + e at step 0, ending with time step 700). The process converges to m at step707.containing deg ( v ) − In order to measure this phenomenon, we introduce the roundness as follows. First, wepartition a configuration into two parts: the outside part with all tiles having deg ( v ) − G = ( V, E ) withsink s and a configuration c : ˜ V → N , let the maximum stable components , MSC ( c ) = { V , V , . . . , V k } , be the connected components of tiles (from ˜ V ) having deg ( v ) − outer tiles is the set outer ( c ) = (cid:91) V i ∈ MSC ( c ) ∃ v ∈ V i : { v,s }∈ E V i and the inner tiles is the set inner ( c ) = ˜ V \ outer ( c ). From this partition of the set oftiles, we are interested in the frontier between the outer and inner tiles, and how closeit is from a perfect circle. This has to do with the coordinates of tiles in the Euclideanspace R , so let us denote coord ( v ) the set of coordinates of the bounds of some tile v ∈ ˜ V , and body ( v ) the subset of R covered by the tile. Our convention regarding thesink s is discussed below. Regarding the circle, let us denote B ( r ) the inside of the circle14igure 14: Stabilization of ( m + e ) ◦ = m on 5 iterations of the substitution from a P2Sun. From left to right, on top are displayed time steps 0 , , ,
150 and at the bottomtime steps 165 , , , m at step 220.of radius r ∈ R centered at coordinate (0 , i.e. B ( r ) = { ( x, y ) ∈ R | (cid:112) x + y ≤ r } . We may therefore denote body ( v ) ∩ B ( r ) = ∅ to state that tile v is entirely outside thecircle of radius r centered at (0 , body ( v ) ⊆ B ( r ) (or equivalently coord ( v ) ⊆ B ( r )since we deal with polygons and the circle is convex) to state that tile v is entirely insidethe same circle. We define the outer radius as the maximum scalar r such that all outertiles are outside the circle of radius r around the origin (0 ,
0) of the tiling, r ( c ) = max { r ∈ R + | ∀ v ∈ outer ( c ) : body ( v ) ∩ B ( r ) = ∅} and the inner radius as the minimum scalar r such that all inner tiles are inside thecircle of radius r around the origin (0 ,
0) of the tiling, r ( c ) = min { r ∈ R + | ∀ v ∈ inner ( c ) : body ( v ) ⊆ B ( r ) } . In order to deal with the case outer ( c ) = ∅ , which may for example be the case on someconfigurations ( m + e ), we add the convention that the sink s is an infinite tile coveringall the space outside the tiling, whose coordinates are the union of all bounds from tileedges adjacent to the sink, and that it always belongs to outer ( c ). As a consequencethe outer radius is upper bounded by the radius of the inscribed circle (inside the finitetiling) with center at (0 , inner ( c ) = ∅ is not problematic since the minimumdefining r ( c ) is taken on R + , and would therefore equal 0 as expected. The inner radiusis upper bounded by the radius of the circumscribed circle (outside the finite tiling) with15igure 15: Examples of roundness measures, circles of outer radius r ( c ) in blue, innerradius r ( c ) in red. Left: square grid of side length 50, after 500 steps from ( m + e )(converges to m in 707 steps), r ( c ) ≈ . − .
224 = 1 . m + e ) (converges to m in 220 steps), r ( c ) ≈ . − .
729 = 1 . , roundness of a configuration c is then measured as the differencebetween the inner and outer radii, r ( c ) = r ( c ) − r ( c ) . Note that 0 ≤ r ( c ), and we have r ( c ) = 0 when the frontier between inner and outertiles is a perfect circle. Two examples of roundness are given on Figure 15. Remark that since all our tiles are polygonal (with three to six sides), we cannot expectto reach roundness 0 (except when the stabilization process has converged to m ). To getsome easy to interpret base values, we consider the diameter of each tile, as the diameterof the circumscribed circle around the tile (smallest radius of a circle having the tileentirely is its interior). See Figure 16. The greatest diameter of some tiling’s tiles can beinterpreted 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), thenall tiles having the center of their circumscribed circle on or outside the circle C can beset as outer tiles, and all tiles having the center of their circumscribed circle inside thecircle C can be set as inner tiles, which results in a roundness smaller than twice theradius of the greatest tile’s radius, i.e. smaller than the greatest tile’s diameter. Weobtain the base roundnesses presented on Table 2.16 √ φ π ) Figure 16: Diameters of all tiles, as the diameter of a circumscribed circle (in purple).All tiles at the same scale.Tiling Squaregrids Triangulargrids Hexagonalgrids P2 tilings P3 tilingsBase roundness √ ≈ .
414 2 √ ≈ . π ) ≈ . Table 2: Base roundnesses as the greatest diameter of some tiling’s tiles.
We now present plots of roundness measured during the stabilization process ( m + e ) ◦ = m on the different tilings considered in this article. We decompose the stabilizationprocess from m + e to m into two phases: ◦ phase 1: the dynamics is erratic, ◦ phase 2: the set of inner tiles slowly shrinks, until reaching inner ( m ) = ∅ .The beginning of phase 2 is defined as the first step such that the inner radius is smalleror equal to the inscribed radius of the tiling (maximum radius of a circle entirely insidethe 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 anyinner tile on the border would otherwise lead to the inner radius being greater than theinscribed radius of the tiling. An important observation is that, in all the experimentspresented in this article and performed during its preparation, once in phase 2 with allborder tiles as outer tiles, then all border tiles remain outer tiles , and that the innerradius remains smaller or equal to the inscribed radius. Two full examples of roundnessplots are given on Figures 17, 19, and their companion Figures 18, 20.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 haveperformed in preparing this article, we have observed similar behaviors for all other sizes. They all remain stable with deg ( v ) − m . Observe that any outer tile receivingsome grain would topple, and that toppling any outer tile would result in toppling the whole maximumstable component it belongs to. quare grid of side length
50 steps0 707units036base 424phase 1 phase 2
Figure 17: Plot of the roundness during the stabilization process ( m + e ) ◦ = m , onthe square grid of side length 50. At each step we plot the inner radius r ( c ) in red,outer radius r ( c ) in blue, and roundness r ( c ) in black. Grid has one row per unit andone column per 10 time steps. For example, at step 3, with m + e ⇒ c ⇒ c ⇒ c ,we observe that r ( c ) = 25 √ ≈ .
355 (the radius of the circumscribed circle aroundthe whole tiling) and that r ( c ) = 25 (the radius of the inscribed circle inside the wholetiling), so that r ( c ) = 25( √ − ≈ .
355 (at step 0 some outer tiles near the x = 0and 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 √
2. See also thecompanion Figure 18.Figure 18: Configurations at steps 0 (this is m + e ), 423 and 424 during the stabilizationprocess ( m + e ) ◦ = m , on the square grid of side length 50. Inner radii in red, outer radiiin blue. One can observe the phase transition occurring at step 424: the inner radiusbecomes smaller or equal to the inscribed radius of the tiling.18 iterations steps0 220units023base 146phase 1 phase 2 Figure 19: Plot of the roundness during the stabilization process ( m + e ) ◦ = m , on thetiling obtained after 5 iterations of the substitution from a P2 Sun. At each step weplot the inner radius r ( c ) in red, outer radius r ( c ) in blue, and roundness r ( c ) in black.Grid has one row per unit and one column per 10 time steps. For example, at step 0we observe that r ( m + e ) ≈ .
940 (the radius of the circumscribed circle around thewhole tiling) and that r ( m + e ) ≈ .
069 (the radius of the inscribed circle inside thewhole tiling), so that r ( m + e ) ≈ . m is reached at step 220. The base roundness of this tiling (dashed) is 2. See also thecompanion Figure 20.Figure 20: Configurations at steps 0 (this is m + e ), 145 and 146 during the stabilizationprocess ( m + e ) ◦ = m , on the tiling obtained after 5 iterations of the substitution froma P2 Sun. Inner radii in red, outer radii in blue. One can observe the phase transitionoccurring at step 146: the inner radius becomes smaller or equal to the inscribed radiusof the tiling. 19 quare grid of side length
500 steps40413 70302units06base
Triangular grid of side length
500 steps50163 69445units08base
Hexagonal grid of side length
200 steps13252 24968units016base
Figure 21: Plots of the roundness during phase 2 of the stabilization process ( m + e ) ◦ = m , on the square grid of side length 500, triangular grid of side length 500, and hexagonalgrid 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 roundnessvalues r ( c ) observed among these 100 steps.20 iterations steps42724 67326units05base P2 Star substitution iterations steps27337 41918units05base P3 Sun substitution iterations steps42067 62333units04base P3 cut and project tiling of size
100 steps39210 66495units04base
Figure 22: Plots of the roundness during phase 2 of the stabilization process ( m + e ) ◦ = m , on the P2 Sun after 11 iterations of the substitution, P2 Star after 11 iterations ofthe substitution, P3 Sun after 10 iterations of the substitution, P3 cut and project tilingof size 100. Grids have one row per unit and one column per 1000 time steps. Twopoints are drawn every 100 time steps: the maximum and minimum roundness values r ( c ) observed among these 100 steps. 21iling Size/iter. Tiles Identity(seconds) Roundness(seconds) Phase 2begin step StabilizationstepSquaregrid 500 250000 825 1294 40413 70302Triangulargrid 500 250000 1143 689 50163 69445Hexagonalgrid 200 120601 271 405 13252 24968P2 Sunsubst. 11 327750 2931 1850 42724 67326P2 Starsubst. 11 234410 944 640 27337 41918P3 Sunsubst. 10 266860 2227 1570 42067 62333P3 cut& project 100 249610 3145 2386 39210 66495Table 3: Data regarding the computations generating the graphics presented on Fig-ures 21 and 22. Running times are given for the computation of identities (from For-mula 1) and roundness measures. The graphics were computed on our personal machines (standard laptops), with simu-lation times and data presented on Table 3. Details on the implementation (no paral-lelization) can be found at https://github.com/huacayacauh/JS-Sandpile/wiki/Roundness . Well, it appears clearly that the circular shapes observed during the stabilization process( m + e ) ◦ = m are not asymptotically approaching perfect circles. Indeed, in all exper-iments conducted on large tilings (Figures 21 and 22), at the beginning of phase 2 theroundness is above the base roundness, whereas the base roundness is an upper boundon the best achievable roundness (see Subsection 4.2).As the circular shapes shrink, it is normal to see the roundness decrease until thevalue 0 on configuration m (all plots reach the minimum roundness 0 at stabilizationstep). 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 phase2 are less round. This is illustrated on Figure 23.Although they do not tend to perfect circles, we may admit that these shapes are closeto be, in regard of tiling sizes (Table 3). Increases of roundness are visible to the nakedeye on hexagonal and triangular grids which deviate largely from a perfect circle, but arehardly noticeable without measurements on other tilings. Figure 24 illustrates this withthe frontier between outer and inner tiles, on the configuration obtained at the beginning22 ide lengthunits06base 50 100 150 200 250 300 350 400 450 500 Figure 23: Maximum roundness r ( c ) measured during phase 2 of the stabilizationprocess ( m + e ) ◦ = m on square grids of various side lengths, illustrating the correlationbetween roundness and size.of phase 2 during the experiments from Subsection 4.3 (since these configurations arequite large, we picture only the frontier in order to reduce the numerical weight of thepresent document). If they are not perfect circles, then we may naturally ask: whatcharacterize these frontier shapes for each tiling?Finally, we see on Figure 24 (first line) that the frontier on grids (respectivelysquare, 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 thenumber of parts where the frontier deviates from a circle (the number of times it goesfrom the inner radius to the outer radius). Theses symmetries are expected, they comefrom the symmetry of the grids and therefore the symmetries of the dynamics. A naturalattempt would be to experiment the roundness of a square grid cropped to a circle, inorder to remove the effect of the border’s shape. Despite the fact that the differencebetween the circumscribed and inscribed radius is smaller than the base roundness, thisdoes not lead to significantly smaller roundness measurements (frontiers are not closerto perfect circles) during phase 2 of the stabilization process ( m + e ) ◦ = m , as shown onFigure 25. It feels that the frontier reflects the anisotropy of the square grid itself ratherthan that of its border’s shape. The experiments presented in this article are reproducible with
JS-Sandpile (links inPreamble). The software implements no parallelization mechanism, which would allowto perform larger simulations ( e.g. using GPU). Nevertheless we believe that this wouldnot lead to qualitatively different observations.We have presented some identity elements on Penrose tilings, revealing no obviousstructure related to these famous aperiodic tilings. Identities are highly sensitive to theshape of the tiling, and it may be the case that other finite croppings of (infinite) Penrosetilings lead to different observations. We tried to build the most “natural” finite Penrosetilings: Suns and Stars obtained by substitution, along with cut and project from a 5-23igure 24: Frontier between inner and outer tiles on the configuration obtained at thebeginning 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 step40413 ( r ( c ) ≈ . r ( c ) ≈ . r ( c ) ≈ . r ( c ) ≈ . r ( c ) ≈ . r ( c ) ≈ . r ( c ) ≈ . adiusunits06base 50 100 150 200 250 Figure 25: Left: maximum roundness r ( c ) measured during phase 2 of the stabilizationprocess ( m + e ) ◦ = m on square grids cropped to circles of various radii. Right: frontierbetween inner and outer tiles on the configuration obtained at the beginning of phase 2(step 32131, r ( c ) ≈ . ), Cayley graphs [53], or hyperbolic planes ( e.g. with Poincar´edisk models).Finally, could the stabilization process ( m + e ) ◦ = m shed some light on the enig-matic 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 breakingsome scientific locks in the domain of sandpiles. The difficulty may be to find constructions from non Wang tiles, because Wang tiles would lead tosquare grids for the sandpile model to play on (as we remove tile decorations). cknowledgments The authors are thankful to Valentin Darrigo for his contributions to
JS-Sandpile , toVictor Poupet for stressing that the apparent isotropy of the ( m + e ) ◦ process is sur-prising at the occasion of a talk given by KP during AUTOMATA’2014 in Himeji, toThomas Fernique for sharing his expertise (and code!) regarding the cut and projectmethod, 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´e,doing an internship at the LIS laboratory (UMR 7020), both in Marseille, France. Thework of KP was funded mainly by his salary as a French State agent and therefore byFrench taxpayers’ taxes, affiliated to Aix-Marseille Univ., Univ. de Toulon, CNRS, LIS,UMR 7020, Marseille, France and Univ. Cˆote d’Azur, CNRS, I3S, UMR 7271, SophiaAntipolis, France. Secondary financial support came from ANR-18-CE40-0002 FANsproject, ECOS-Sud C16E01 project, and STIC AmSud CoDANet 19-STIC-03 (CampusFrance 43478PD) project. References [1] M. Baake, M. Scholottmann, and P. D. Jarvis. Quasiperiodic tilings with tenfoldsymmetry 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 ofthe 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 ofTilings. In
Journ´ees Automates Cellulaires 2010 , pages 54–64, 2010.[5] R. Berger. The undecidability of the domino problem.
Memoirs of the AmericanMathematical Society , 66, 1966.[6] H. Cairns. Some halting problems for abelian sandpiles are undecidable in dimensionthree.
SIAM Journal on Discrete Mathematics , 32(4):2636–2666, 2018.[7] J. Cervelle and B. Durand. Tilings: recursivity and regularity.
Theoretical ComputerScience , 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 2Dcellular automata.
Theoretical Computer Science , 218(2):347–417, 1999.2610] D. Dhar. Self-organized critical state of sandpile automaton models.
Physical ReviewLetters , 64:1613–1616, 1990.[11] D. Dhar, P. Ruelle, S. Sen, and D. N. Verma. Algebraic aspects of abelian sandpilemodels.
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.
TheoreticalComputer Science , 205(1-2):183–193, 1998.[14] V. G. Fast and I. R. Efimov. Stability of vortex rotation in an excitable cellularmedium.
Physica D: Nonlinear Phenomena , 49(1):75–81, 1991.[15] T. Fernique.
Pavages, Fractions Continues et G´eom´etrie Discr`ete . PhD thesis,Universit´e Montpellier II, 2007.[16] E. Formenti, E. Goles, and B. Martin. Computational complexity of avalanches inthe Kadanoff sandpile model.
Fundamenta Informaticae , 115(1):107–124, 2012.[17] E. Formenti, B. Masson, and T. Pisokas. Advances in symmetric sandpiles.
Funda-menta Informaticae , 76(1-2):91–112, 2007.[18] E. Formenti, K. Perrot, and E. R´emila. Computational complexity of theavalanche problem on one dimensional Kadanoff sandpiles. In
Proceedings of AU-TOMATA’2014 , volume 8996 of
LNCS , pages 21–30, 2014.[19] E. Formenti, K. Perrot, and E. R´emila. Computational complexity of the avalancheproblem 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 theparallel symmetric sandpile model.
Theoretical Computer Science , 533:1–14, 2014.[21] Enrico Formenti and K´evin 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.
The-oretical Computer Science , 369(1-3):463–469, 2006.[23] E. Goles. Sand pile automata.
Annales de l’institut Henri Poincar´e (A) Physiqueth´eorique , 56(1):75–90, 1992.[24] E. Goles and M. Kiwi. Games on line graphs and sand piles.
Theoretical ComputerScience , 115(2):321–349, 1993.[25] E. Goles, M. Latapy, C. Magnien, M. Morvan, and H. D. Phan. Sandpile modelsand lattices: a comprehensive survey.
Theoretical Computer Science , 322(2):383–407, 2004. 2726] E. Goles, D. Maldonado, P. Montealegre, and N. Ollinger. On the computationalcomplexity of the freezing non-strict majority automata. In
Proceedings of AU-TOMATA’2017 , pages 109–119, 2017.[27] E. Goles and M. Margenstern. Universality of the chip-firing game.
TheoreticalComputer Science , 172(1-2):121–134, 1997.[28] E. Goles and P. Montealegre. Computational complexity of threshold automatanetworks under different updating schemes.
Theoretical Computer Science , 559:3–19, 2014.[29] E. Goles and P. Montealegre. A fast parallel algorithm for the robust predictionof 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 SystemSciences , 91:1–32, 2017.[31] E. Goles, P. Montealegre-Barba, and I. Todinca. The complexity of the bootstrapingpercolation and other problems.
Theoretical Computer Science , 504:73–82, 2013.[32] E. Goles, M. Morvan, and H. D. Phan. Sandpiles and order structure of integerpartitions.
Discrete Applied Mathematics , 117(1–3):51–64, 2002.[33] B. Gr¨unbaum 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 inavalanches.
Physical Review A , 39(12):6524–6537, 1989.[35] L. Levine, W. Pegden, and C. K. Smart. Apollonian structure in the abelian sand-pile.
Geometric and Functional Analysis , 26:306–336, 2016.[36] L. Levine, W. Pegden, and C. K. Smart. The apollonian structure of integer super-harmonic matrices.
Annals of Mathematics , 186(1):1–67, 2017.[37] L. Levine and Y. Peres. Laplacian growth, sandpiles, and scaling limits.
Bulletinof the American Mathematical Society , 54(3):355–382, 2017.[38] M. Marek. Grid anisotropy reduction for simulation of growth processes with cellularautomaton.
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.2841] C. Moore. Majority-vote cellular automata, Ising dynamics, and P-completeness.[42] C. Moore and M. Nilsson. The computational complexity of sandpiles.
Journal ofStatistical 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 , volume10875 of
LNCS , pages 127–142, 2018.[44] W. Pegden and C. K. Smart. Convergence of the abelian sandpile.
Duke Mathe-matical Journal , 162(4):627–642, 03 2013.[45] W. Pegden and C. K. Smart. Stability of patterns in the abelian sandpile.
AnnalesHenri Poincar´e , 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 Mathe-matical Intelligencer 2 , pages 32—-37, 1979.[48] K. Perrot.
Les piles de sable Kadanoff . PhD thesis, ´Ecole normale sup´erieure deLyon, 2013.[49] K. Perrot, H. D. Phan, and V. T. Pham. On the set of Fixed Points of the ParallelSymmetric Sand Pile Model. In
Proceedings AUTOMATA’2011 , DMTCS, pages17–28. Open Publishing Association, 2011.[50] H. D. Phan.
Structures ordonn´ees et dynamiques de piles de sable . PhD thesis,Universit´e 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.
Inven-tiones mathematicae , 12:177–209, 1971.[53] Zs. Roka.
Automates cellulaires sur graphes de Cayley . PhD thesis, ´Ecole NormaleSup´erieure de Lyon, 1994.[54] H. E. Schepers and M. Markus. Two types of performance of an isotropic cellularautomaton: stationary (Turing) patterns and spiral waves.
Physica A: StatisticalMechanics and its Applications , 188(1):337–343, 1992.[55] B. Sch¨onfisch. Anisotropy in cellular automata.
Biosystems , 41(1):29–41, 1997.[56] G. Ch. Sirakoulis, I. Karafyllidis, and A. Thanailakis. A cellular automaton forthe propagation of circular fronts and its applications.
Engineering Applications ofArtificial Intelligence , 18(6):731–744, 2005.2957] H. Wang. Proving theorems by pattern recognition — II.
The Bell System TechnicalJournal , 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 incellular automaton models of excitable media.