Multigrid Generation of Aperiodic Tilings

Basic Concepts

The Penrose tiling is the classic aperiodic tiling of the plane. In its simplest form, it consists of two tile shapes - a thick rhombus with angles of 72 and 108 degrees and a thin rhombus with 36 and 144 degree angles. Both rhombs have the same length sides. Penrose marked the rhombs so that they only fit together in certain orientations, but there are are also complex non-local requirements which make fitting together a large number of tiles a rather hit-or-miss affair. However, de Bruijn discovered a very clever deterministic way to fit together an arbitrarily large number of tiles without ever making any mistakes. Here, we will explore generalizations of de Bruijn's pentagrid idea, leading to an algorithm to draw Penrose and other similar tilings.

line intersections generate dual rhomb tiling

De Bruijn noticed that any set of intersecting lines in the plane defines a dual tiling as follows: Imagine an infinitessimal rhombus with its centroid at each intersection point and its edges perpendicular to the lines. If you inflate all of these to have a common side length, then you can use the sequence of intersection points along any one line to tell you the order to fit the corresponding rhomb edges together to completely cover a connected convex region of the plane. (If more than two lines intersect at a point, you can use the equilateral polygon whose edges are perpendicular to the lines as a single tile corresponding to that point - or you can decompose that even-sided convex polygon into rhombs by infinitessimally displacing the lines.)

(See adamponting.com/de-bruijn-tilings for nicer examples.)

To recap: Any set of intersecting lines determines a dual rhomb tiling. Each rhomb corresponds to one intersection point, and each vertex in the tiling corresponds to a polygonal region bounded by the lines. All the rhomb tile edges have equal lengths and the edges of each rhomb are perpendicular to the two lines of its corresponding intersection point. Walking from point to point around a polygonal region in the plane of the lines corresponds to walking from rhomb to rhomb around a vertex in the tiling. And walking around the four regions surrounding each intersection point corresponds to walking around the four corners of one rhomb in the tiling.

Now de Bruijn looked specifically at five infinite grids of equally spaced parallel lines, with each set of grid lines parallel to one edge of a regular pentagon, and discovered that not only did the dual rhomb tilings satisfy the Penrose matching conditions, but also that every infinite Penrose tiling was dual to such a “pentagrid”. Here we consider $n$ sets of parallel grid lines, first with arbitrary directions and line spacings, then with equal angles between the grids and equal line spacings more analogous to pentagrids.

Generic Algorithm for Finite Number of Lines

Before we move on to infinite sets of parallel grid lines, consider a collection of $n$ random lines $\ZETAH_j\cdot\VEC{x}+\gamma_j=0.$ Generally we have $n(n-1)/2$ intersection points - one intersection for each pair of lines. An algorithm to draw the dual rhomb tiling is fairly simple:

  1. Loop over the $n$ lines, with line index $j$.
  2. Find the $n-1$ intersections with the other lines, with line index $k\ne j$.
  3. Each intersection $(j,k)$ with $j<k$ will correspond to a rhomb with corners $\VEC{z}_{jk},$ $\VEC{z}_{jk}+\ZETAH_j,$ $\VEC{z}_{jk}+\ZETAH_j+\ZETAH_k,$ and $\VEC{z}_{jk}+\ZETAH_k.$
  4. Sort the intersections with $j$ in order of $\ZETAH_j\cdot\VEC{x}.$
  5. If $j=0$, find the intersection $k$ in the middle of the sorted list and draw that rhomb at $\VEC{z}_{jk}=(0,0).$ For $j>0,$ the process is the same, except you begin from the $k=0$ intersection.
  6. Step forward in the sorted list to $k',$ and if $k'>j,$ draw that rhomb at $\VEC{z}_{jk'}=\VEC{z}_{jk}+\ZETAH_k,$ otherwise just move to that new position.
  7. After you reach the end of the sorted list, return to the initial $k$ and step backward, this time drawing or skipping the rhombs at $\VEC{z}_{jk'}=\VEC{z}_{jk}-\ZETAH_{k'}.$
  8. Note that for $j=0$ you must keep the list of rhomb positions $\VEC{z}_{0k}$, indexed by $k.$ You can forget the rhomb locations for $j>0,$ just keeping track of the current $\VEC{z}_{jk}.$

The more structured infinite grids we consider next eliminate the need for the sorting. However, the sort operations do not make these finite tilings particularly more expensive to draw than the grids we turn to now.

Arbitrary Grid Directions and Spacings

Consider an arbitrary set of non-zero, non-parallel two dimensional vectors $\{\ZETA_j\}$ for $j = 1, 2, ..., n,$ and $n$ grids of equally spaced parallel lines \[ \ZETA_j\cdot\VEC{x} + \gamma_j = k_j, \] where $k_j$ runs through all integer values and $\gamma_j$ is an arbitrary constant determining the phase of the grid lines relative to the origin. Note that changing $\gamma_j$ by any integer value merely relabels the $k_j$ for all of the grid lines without changing the grid itself. Any pair $r, s$ of grids generates a regular lattice of intersection points $\VEC{x}(k_r, k_s)$ labeled by the line indices $(k_r, k_s).$

Although each $(r, s)$ pair produces a regular set of lattice points, when the $\ZETA_j$ are incommensurate the lattices corresponding to different pairs never synchronize, so the full pattern of $n(n-1)/2$ lattices is aperiodic. Nevertheless, we can construct the dual tiling explicitly.

Every point $\VEC{x}$ in the multigrid plane falls between one pair of $j$-lines, $k_j=K_j(\VEC{x})-1$ and $k_j=K_j(\VEC{x})$, which we can write using ceiling brackets: \[ K_j(\VEC{x}) = \lceil \ZETA_j\cdot\VEC{x} + \gamma_j \rceil. \] We also define the fractional distance from the point $\VEC{x}$ to the line $K_j(\VEC{x})$, $0 <= \lambda_j < 1,$ where \[ \lambda_j(\VEC{x}) = K_j(\VEC{x}) - \ZETA_j\cdot\VEC{x} - \gamma_j. \]

Next, define the function \[ f(\VEC{x}) = \sum_j K_j(\VEC{x})\ZETAH_j, \] where $\ZETAH_j$ is the unit vector in the direction of $\ZETA_j.$ The function $f$, like the $n$ functions $K_j$, is constant within every polygonal region of the multigrid. Therefore, the value of $f$ in each region corresponds to a vertex in the dual tiling. For vertex $(k_r, k_s)$ in the $rs$ lattice, call the value of $f$ \[ \VEC{z}(k_r, k_s) = \sum_j K_j(\VEC{x}(k_r,k_s))\ZETAH_j. \] Furthermore, as you walk around the intersection point $\VEC{x}(k_r,k_s),$ $f(\VEC{x})$ successively takes the four values \[ \VEC{z}(k_r, k_s),\quad \VEC{z}(k_r+1, k_s),\quad \VEC{z}(k_r+1, k_s+1),\quad \VEC{z}(k_r, k_s+1). \] Since crossing an $r$ (or $s$) line in the multigrid changes just $K_r$ (or $K_s$) by $\pm 1,$ these four values of $f$ are in fact the four corners of a unit-sided rhombus: \[ \VEC{z}(k_r, k_s),\quad \VEC{z}(k_r, k_s)+\ZETAH_r,\quad \VEC{z}(k_r, k_s)+\ZETAH_r+\ZETAH_s,\quad \VEC{z}(k_r, k_s)+\ZETAH_s \] Therefore, not only do the values of $f(\VEC{x})$ correspond to the vertices in a rhomb tiling, we can identify them as the vertices in a rhomb tiling dual to the multigrid! All the rhomb sides in this tiling have unit length, and the number of different rhomb shapes corresponds to the number of different angles between $\ZETAH_r$ and $\ZETAH_s$ for all $rs$ pairs.

Furthermore, we have discovered an algorithm to actually draw the dual tiling given the multigrid: Loop over all $n(n-1)/2$ of the $rs$ pairs, and draw the rhombs for each $(k_r, k_s)$ pair where that intersection point lies in some bounded region of the multigrid space.

Now rewrite the formula for $\VEC{z}(k_r, k_s)$ (the corner of the rhomb associated with intersection point $k_r, k_s$ that has the smallest values of $k_r$ and $k_s$) in terms of the remainders $\lambda_j(k_r, k_s) = \lambda_j(\VEC{x}(k_r, k_s))$. By definition $K_j(\VEC{x}(k_r, k_s)) = \ZETA_j\cdot\VEC{x}(k_r,k_s) + \gamma_j + \lambda_j(k_r,k_s),$ so \[ \VEC{z}(k_r, k_s) = \VEC{M}\cdot\VEC{x}(k_r,k_s) + \ETA + \sum_j\lambda_j(k_r,k_s)\ZETAH_j \] where we have defined the symmetric $2\times 2$ matrix $\VEC{M}$ as the sum of outer products \[ \VEC{M} = \sum_j\ZETA_j\,\ZETAH_j, \] and the vector $\ETA$ as \[ \ETA = \sum_j\gamma_j\,\ZETAH_j. \]

That is, the point $\VEC{z}$ in the tiling plane is always close to the point $\VEC{M}\cdot\VEC{x}(k_r,k_s) + \ETA + \sum_j\ZETAH_j/2$, where we have added the final term as the mean of the possible range of values of the $\lambda_j$ sum. Since we are primarily interested in the rhomb tiling, we can consider the sets of grid lines with the linear transform $\VEC{M}$ applied. This transformed multigrid will still have the same topology of intersection points as the original grid. However, instead of a rhomb associated with each intersection point, in general the transform $\VEC{M}$ will distort it into a parallelogram. In terms of $\VEC{x}'=\VEC{M}\cdot\VEC{x}$ the equations for the grid lines are \[ \PSI_j\cdot\VEC{x}' + \gamma_j = k_j \] where $\PSI_j = \VEC{M}^{-1}\cdot\ZETA_j.$ We also translate the definition of $\VEC{z}(k_r,k_s)$ to make it line up with $\VEC{x}'$ \[ \VEC{z}(k_r, k_s) = \sum_j K_j(\VEC{x}(k_r,k_s))\ZETAH_j - \ETA - \tfrac{1}{2}\sum_j\ZETAH_j, \] so that \[ \VEC{z}(k_r, k_s) = \VEC{x}'(k_r,k_s) + \sum_j(\lambda_j(k_r,k_s)-\tfrac{1}{2})\ZETAH_j. \]

Since the $\lambda_j$ are always between 0 and 1, we see that the center of the $(k_r,k_s)$ rhomb is always near its dual point $\VEC{x}'(k_r,k_s).$ Furthermore, $\lambda_r(k_r,k_s) =\lambda_s(k_r,k_s)=0,$ so there are only $n-2$ non-zero terms in the displacement $\sum_j\lambda_j(k_r,k_s)\ZETAH_j.$

Note that the grid intersection points can be written explicitly as \[ \VEC{x}'(k_r, k_s) = (k_s - \gamma_s)\PHI_{rs} + (k_r - \gamma_r)\PHI_{sr} \] where $\PHI_{rs} = \PSI_{r\perp} / \PSI_r \times \PSI_s,$ and $\PSI_{r\perp}$ is $\PSI_r$ rotated ninety degrees counterclockwise so that $\PSI_{r\perp} \cdot \PSI_s = \PSI_r \times \PSI_s.$ Beware that a set of $k_j$ values uniquely identifies only a zone in the multigrid (or a vertex in the tiling) - not an intersection point in the multigrid. A set of $k_j$ values may be shared by different intersection points in separate $(r, s)$ lattices, or by the corresponding different rhombs in the tiling.

Since each rhomb shape corresponds to an $(r, s)$ pair, the density of intersection points of that lattice determines the density of that rhomb shape. The ratio of these densities, in turn, must equal the ratio of numbers of the corresponding rhomb shapes. The irrationality of these density ratios guarantees a tiling is aperiodic.

Equally Spaced Multigrids

When the angles are equally spaced - at $n$-th roots of unity when $n$ is odd, or $2n$-th roots of unity when $n$ is even - and all $n$ grid spacings are equal, the only real simplification in this picture is that the matrix $\VEC{M}$ is proportional to the identity matrix. In fact, $\VEC{M}=\tfrac{n}{2}\VEC{I}$ so that $\PSI_j=\tfrac{2}{n}\ZETAH_j.$ (We assume $\ZETA_j=\ZETAH_j$ here.) This means that $\VEC{x}'=\tfrac{n}{2}\VEC{x}$ is an undistorted version of the original multigrid.

For the case of odd $n$, something else interesting happens. In this case, $\sum_j\ZETA_j=0.$ If we additionally restrict the $\gamma_j$ to have the property $\sum_j\gamma_j=0,$ then since $K_j(\VEC{x})$ is an integer, \[ \begin{eqnarray*} \sum_j K_j(\VEC{x}) &=& \sum_j(\ZETA_j\cdot\VEC{x} + \gamma_j) + \sum_j\lambda_j(\VEC{x}) \\ &=& \sum_j\lambda_j(\VEC{x}) \\ &=& I(\VEC{x}) \end{eqnarray*} \] is also an integer for every point $\VEC{x}$ in the multigrid plane. Since $0\le\lambda_j\lt 1,$ we have $I(\VEC{x})\lt n.$ Furthermore, assuming all $n$ lines never meet at any point, $I(\VEC{x})\gt 0.$ The function $I(\VEC{x})$ is called the index of a point.

For the $n=5$ case de Bruijn first studied, this $\sum_j\gamma_j=0$ restriction (equivalent to summing to any integer value, since shifting the $k_j$ line labels does not change the grid) plays an important role in his proof that the pentagrid construction is identical to Penrose's original edge matching construction. This proof is based on the integer values of the index, which he can use to label the vertices in the tiling. Those vertex labels can be used to reconstruct the edge matching conditions. If $\sum_j\gamma_j$ has some non-zero, non-integer value, it is unclear what happens - the tilings derived from such a pentagrid still look like Penrose tilings, but perhaps the edge matching conditions can occasionally fail.

What is the analogous vertex labeling function for even $n$? When you impose the $\sum_j\gamma_j=0$ restriction for even $n$, the index $I(\VEC{x})$ is again between $0$ and $n$ (non-inclusive), but it is no longer restricted to an integer value, since $\sum_j\ZETA_j\ne 0$ in this case. However, the values of $I(\VEC{x})$ at the four corners of every rhomb still differ by exactly $1$ for corners connected by a rhomb edge.

The quantities $\sum_j\gamma_j$ and $\ETA=\sum_j\gamma_j\,\ZETA_j$ are the zeroth and first Fourier components (interpreted as complex numbers) of the real sequence $\gamma_j.$ Specifically, this is a discrete Fourier transform on $n$ elements for odd $n$, or $2n$ elements for even $n$. In the even $n$ case, the grid line equations $\ZETA_j\cdot\VEC{x} +\gamma_j=k_j$ show that the second half of the $2n$ DFT components must be $\gamma_{n+j}=-\gamma_j,$ since $\ZETA_{n+j}=-\ZETA_j.$

The $k$th element of the DFT is \[ \XI_k = \sum_{j=0}^{N-1}\gamma_j\ZETA_{kj}, \] where $n=n$ for odd $n$, $N=2N$ for even $n.$ Since the $\gamma_j$ are real, $\XI_{-k}=\overline{\XI_{k}}.$ (Here the overline means complex conjugate.) For even $n,$ $\XI_k=$ for all even $k$ since $\gamma_{n+j}=-\gamma_j.$ Furthermore, for odd $k$, the second set of $n$ terms in the sum on $j$ exactly equals the first $n$ terms, because $\ZETA_n=-1$ for even $n.$ Therefore, we will renormalize the DFT coefficients for even $n$ by a factor of two so that for all $n$ we define \[ \XI_k = \sum_{j=0}^{n-1}\gamma_j\ZETA_{kj}. \] For even $n,$ this formula is only valid for odd $k$; for even $k,$ $\XI_k=0$ from our consideration of the full DFT sum. For odd $n$ all values of $k$ for $\XI_k$ are valid, but we need only consider $k\lt n/2$ because the other half are simply complex conjugates.

We have already noted that $\ETA=\XI_1$ represents a simple displacement of the grid, while the condition that $\XI_0=0$ was important in our discussion of the index function.

The inverse DFT for odd $n$ is \[ \begin{eqnarray*} \gamma_j &=& \tfrac{1}{n}\sum_{k=0}^{(n-1)/2}(\XI_k\ZETA_{-kj} +\overline{\XI_k}\ZETA_{kj}) \\ &=& \tfrac{1}{n}(\XI_0 + 2\sum_{k=1}^{(n-1)/2}\XI_k\cdot\ZETA_{kj}). \end{eqnarray*} \] Only half of the $\XI_k$ are required in the sum; their real and imaginary parts contain twice as much information as the original $\gamma_k$, and $\XI_0$ is real.

The inverse DFT for even $n$ is (remembering we have reduced the DFT $\XI_k$ by a factor of two, so the denominator is $n$ instead of $N=2n$): \[ \begin{eqnarray*} \gamma_j &=& \tfrac{1}{n}\sum_{k=0}^{n-1}(\XI_k\ZETA_{-kj} +\overline{\XI_k}\ZETA_{kj}) \\ &=& \tfrac{2}{n}\sum_{k=1}^{n-1}\XI_k\cdot\ZETA_{kj}. \end{eqnarray*} \] Remember that $\XI_k=0$ for even $k$ (including $k=0$), so the sum contains only odd $k.$ With their real and imaginary parts these odd $\XI_k$ have the same number of degrees of freedom as the original $\gamma_j.$

Boyle has an extraordinarily detailed analysis of the $n$ equals 4, 5, and 6 cases, in which he proves that the self-similar Ammann-Beenker, Penrose, and Socolar tilings, respectively, can be represented by multigrids precisely when the $\gamma_j$ have $\XI_k$ which are non-zero only for $k=1$ and $k=(n-1)/2$ for odd $n$ ($n=5$) or only for $k=1$ and $k=n-1$ for even $k.$ Here, I adopt this constraint for initializing equally spaced multigrids for all values of $n.$ Adding non-zero $\XI_k$ for other values of $k$ than these often produces distinctly different kinds of vertices in the tiling pattern. Somehow these Fourier decompositions of the $\gamma_j$ seem to be significant beyond the $n=6$ case.

The case of $n$ being a multiple of three is an example of how other Fourier coefficients affect the kinds of vertices which appear in the tilings. These tilings all include a 60 degree rhomb. It is not difficult to prove that as long as $\XI_k=0$ for all $k$ which are multiples of three (again assuming $n$ is a multiple of three), all the 60 degree rhombs in the tiling will meet at one of their 120 degree angles to form hexagons.

The case of $n=2$ is just a periodic square tiling. The case of $n=3$ is always a periodic triangular multigrid producing a periodc tiling of 60 degree rhombs. Apparently, the reason that $n=4$ and $n=6$ have their additional self-similar structure is that they represent an overlay of exactly two of these periodic lattices. The $n=5$ case also has only two different lattice angles. Larger $n$ have more different kinds of periodic lattices overlaid.