Three variations on a theme by Fibonacci

Several variants of the classic Fibonacci inflation tiling are considered in an illustrative fashion, in one and in two dimensions, with an eye on changes in the spectrum. In one dimension, we consider extension mechanisms of deterministic and of stochastic nature, while we look at direct product variations in a planar extension. For the pure point part, we systematically employ a cocycle approach that is based on the underlying renormalisation structure and allows explicit calculations also in cases where one meets regular model sets with Rauzy fractals as windows.


Introduction
The mathematical theory of aperiodic order profits enormously from two construction principles, namely the cut and project method and the substitution or inflation method; see [7,13] and references therein for general background. While the understanding of cut and project sets has reached a rather satisfactory level, this is less so for primitive inflation tilings. In particular, the classification of their spectra is as yet incomplete. While the still open Pisot substitution conjecture [2,21] marks the frontier for one-dimensional systems, the situation is worse in two or more dimensions. Here, it is already less clear what an appropriate conjecture could be, because there are rather intricate additional constraints of geometric origin that do not exist in one dimension.
In this exposition, we reconsider the simplest and best-studied example, the Fibonacci chain, and investigate modifications as well as extensions to planar systems. On the line, we demonstrate two mechanisms that add continuous spectral components due to disorder, either of random or of deterministic type. In the case of planar tilings, we are particularly interested in geometric variations that change the system topologically, but not metrically, which is to say that we are probing the stability of spectral properties. Our starting point is the self-similar Fibonacci tiling dynamical system (TDS) in one dimension, as defined by the primitive inflation rule ̺ : a → ab , b → a, with a and b considered as tiles (or intervals) of length τ = 1 2 1 + which is primitive, with Perron-Frobenius (PF) eigenvalue τ and corresponding left and right eigenvectors, in Dirac notation, u| = τ + 2 5 τ, 1 and |v = τ −1 , τ −2 T .
If Y denotes the (compact) tiling hull defined by ̺, the topological dynamical system (Y, R) is strictly ergodic; compare [11,Ch. 5]. More precisely, Y can be constructed as the orbit closure of a bi-infinite Fibonacci tiling, for which one usually employs one of the two fixed points of ̺ 2 , with core a|a or b|a; see [7,Ex. 4.6] for details and [23] for general background. Since each fixed point is repetitive, the resulting dynamical system is minimal. The unique probability measure on Y is the patch frequency measure, denoted by µ. Via the left endpoints of the intervals, one can consider the elements of Y either as tilings or as Delone sets, which are two viewpoints that we will tacitly identify. This is justified by the fact that the two structures are mutually locally derivable (MLD) from one another; see [7,Sec. 5.2] for background. Note that, in our setting, each such Delone set has density (τ + 2)/5 = τ / √ 5. The topological dynamical system (Y, R) is uniquely ergodic and has pure point spectrum in the measure-theoretic sense, employing the Koopman operator on the Hilbert space L 2 (Y, µ). All eigenfunctions are continuous [15,Thm. 1.4], though the system itself is not equicontinuous. In the length scale chosen, the spectrum is , in addtive notation, as can easily be extracted from the description of Y by the projection method; see [7,Sec. 9.4.1]. The understanding of this example is based on its simultaneous description as an inflation tiling and as a cut and project set. We will recall the details in Section 2, based on the Minkowski embedding of Z[τ ] as a lattice in R 2 .
Below, we go through a number of modifications or variations, starting with one-dimensional examples with some added structure, random or deterministic, to show origins of continuous spectral components. Then, we double the dimension and consider Fibonacci-type tilings of the plane, and some of their rearrangements via modified inflation rules. Here, we present two rather typical phenomena, namely the robustness of pure point spectra in the Pisot case and the emergence of Rauzy fractals already in simple situations.
Our exposition has a somewhat informal character, aiming at an illustration of several possible phenomena via simple examples. To do so, we present the various derivations first, and then summarise our results in formal theorems. We begin with a brief summary of the Fibonacci TDS in Section 2, followed by variations of it that add absolutely continuous or singular continuous components (Section 3). Then, in Section 4, we define the natural direct product Fibonacci TDS in the plane, and derive its properties. Finally, we investigate direct product variations (DPVs) in Section 5, which show robustness of the spectral type.

The Fibonacci chain and its spectral properties
The fact that (Y, R) has pure point (dynamical) spectrum is equivalent to the statement that the Fibonacci chain has pure point diffraction; see [9] and references therein for the general theory. Here, for any Delone set Λ ∈ Y, one defines the measure ω = δ Λ := x∈Λ δ x , called the Dirac comb of Λ. Then, ω leads to the autocorrelation γ = ω ⊛ ω, where ⊛ denotes volume-averaged (or Eberlein) convolution [7,Sec. 8.8], and to the diffraction measure γ, which is the (existing) Fourier transform of the autocorrelation. One obtains the pure point measure where the amplitudes, or Fourier-Bohr (FB) coefficients, are given by The limits exist uniformly, and one has which means that, for any k for which the coefficient is non-trivial, Λ → A Λ (k) defines an eigenfunction of our system, with eigenvalue k in additive notation. Since eigenfunctions for primitive inflations tilings are continuous, their knowledge on the defining orbit suffices to determine them. For later convenience, we consider the two points sets Λ a and Λ b of left endpoints of tiles of types a and b separately, with Λ = Λ a∪ Λ b , and define corresponding amplitudes A Λa and When Λ is one the two fixed points mentioned above and k ∈ L ⊛ , the amplitudes are given [7,Sec. 9.4

.1] by
are the closures of the windows for the point sets Λ a and Λ b in the projection formalism. For the sum of the amplitudes, the formulas simplify to where sinc(x) = sin(x)/x. Consequently, the diffraction intensities are Note that the intensity function, and thus the diffraction measure γ, is the same for all Λ ∈ Y, while the amplitudes will usually differ by a phase factor for distinct elements.
More generally, if we introduce two (in general complex) weights for the two different endpoints (u a and u b say, so ω = u a δ Λa + u b δ Λ b ), we obtain the intensity as Let us sum up these well-known spectral properties [20,7] as follows.
Theorem 2.1. The Fibonacci dynamical system (Y, R), in its geometric realisation as described above, is strictly ergodic and has pure point spectrum, both in the diffraction and in the dynamical sense. For given weights u a , u b ∈ C, the diffraction measure is given by with the intensities I(k) from Eq. (2.4). The Fourier module is L ⊛ = Z[τ ]/ √ 5 and agrees with the dynamical spectrum of (Y, R) in additive notation.
While the relation between the FB coefficients and the Fourier transform of the compact windows holds for regular model sets in general, it is impossible to compute the coefficients this way if the windows are compact sets with fractal boundaries. Let us therefore explain a different approach that will also work in such, more complicated situations.
With σ = τ ⋆ = 1 − τ , the inflation structure induces a relation between the windows that, in terms of their characteristic functions, reads as can easily be verified for the given windows. Observing that 1 σWa∪σW b = 1 σWa + 1 σW b holds as an equation of L 1 -functions, an application of the inverse Fourier transform gives By an elementary calculation, one finds which holds for arbitrary α, β ∈ R with α = 0 and any compact set K ⊂ R.
Defining h a,b =1 W a,b , an application of (2.6) to Eq. (2.5) results in where B is related to the Fourier matrix from the renormalisation approach [5,4] by first taking the ⋆-map of the set-valued displacement matrix T =

{0} {0} {τ } ∅
and then its (inverse) Fourier transform. For this reason, we call it the internal Fourier matrix [8]. The latter can also be viewed as emerging from the commutative diagram where F denotes the Fourier transform of a matrix of (finite) Dirac combs and ⋆ the induced mapping on the level of the Fourier matrices. Using the notation |h = (h a , h b ) T and applying the above iteration n times leads to which can be interpreted as a singular case of a transfer matrix approach to a cocycle; compare [1,18]. In particular, B (1) = B and B (n) (0) = M n for all n ∈ N, where M is the substitution matrix from (1.1). Note that B (n) (y) defines a matrix cocycle, called the internal cocycle, which by (2.8) is related to the usual inflation cocycle by an application of the ⋆-map to the displacement matrices of the powers of the inflation rule [6,8].
It is not difficult to prove that |h(y) = C(y)|h(0) , where exists pointwise for every y ∈ R. In fact, one has compact convergence, which implies that C(y) is continuous [8,Thm. 4.6 and Cor. 4.7]. Clearly, for any m, n ∈ N, one has Employing this relation with m = 1, letting n → ∞, and observing |σ| = τ −1 , one obtains This relation implies that each row of C(y) is a multiple of u|, so we get C(y) = |c(y) u| with |c(y) = c a (y), c b (y) T . Since C(0) = P with the projector P = |v u| from (1.2), we have |c(0) = |v . As |h(y) = |c(y) u|h(0) , where |h(0) = τ |v follows from a simple calculation, we get |h(y) = τ |c(y) , and the inverse Fourier transforms of the windows are encoded in the matrix C. For the Fibonacci case at hand, we can explicitly calculate |c(y) through the Fourier transforms of the known windows to be from which one can explicitly check that there is no y ∈ R for which both functions vanish simultaneously. Consequently, C(y) is always a matrix of rank 1.
In other situations, where no explicit formula for the Fourier transforms is available, one can replace C(y) by |σ| n B (n) (y) or its analogue for a sufficiently large n, subject to the condition that C(y) is approximated sufficiently well and that |σ| n y (or its analogue) is close enough to 0. This works because the windows are compact sets, so that their Fourier transforms are continuous functions. The convergence of this approximation is exponentially fast; see [8] for details and an extension of the cocycle approach to more general inflation systems.
With our previous relation (2.3), for k ∈ L ⊛ , the FB amplitudes are which means that they can now be calculated via C as well. Unless stated otherwise, numerical calculations and illustrations will always be based on the cocycle approach due to its superior speed and accuracy in the presence of complex windows.

First variation -some randomness or disorder
Given a Fibonacci tiling from the hull Y, we now introduce some (uncorrelated) disorder into the a positions, by randomly assigning two different labels or weights. This way, we address and answer a question by Strungaru [22] on the spectral consequences of this type of modification. Consider the situation where tiles of type a are independently replaced by a tile a of the same length with probability q, and kept with probability p = 1 − q. By a simple calculation, whenever p ∈ (0, 1), this leads to a shift space with topological entropy log(2)/τ , calculated per tile (rather than per unit length).
By the strong law of large numbers, using the results of [3], the pair correlation coefficients for this modified system can be expressed in terms of the original coefficients as follows, . Now, we have three different types of points. Using weights u a , u a and u b , we thus obtain the autocorrelation Setting v a = pu a + qu a and v b = u b , an explicit computation shows that, for z = 0, this becomes γ Taking the Fourier transform, we see that the additional point measure at z = 0 gives rise to an absolutely continuous component, and we obtain the following result, which can also be seen as a special case of the random cluster model treated in [3].
Theorem 3.1. If the Fibonacci chain is randomised at the a-positions with a binary Bernoulli process with probabilities p ∈ [0, 1] and q = 1 − p as described above, the diffraction measure is almost surely given by where λ L denotes Lebesgue measure on R. In particular, for pq = 0 and u a = u a , it is a sum of a pure point and an absolutely continuous measure.
Clearly, if u a = u a or pq = 0, the absolutely continuous component vanishes, and we recover the result for the perfect Fibonacci chain.
Inspired by [16,Sec. 8] and [8], let us next extend the Fibonacci substitution to a four-letter substitution with bar-swap symmetry [5], namely The substitution matrix is and the PF eigenvectors u| = 2(3τ +4) 19 (τ, τ, 1, 1) and |v = 2−τ 2 (τ, τ, 1, 1) T . Normalisation is again such that u|v = 1|v = 1. Let us choose natural interval lengths τ for a, a and 1 for b, b, which means that we use the same setting as for the perfect Fibonacci chain, including the ring of integers, Z[τ ], and its Minkowski embedding L from (2.2). The one-sided fixed point equations for the Delone sets of left endpoints read Applying the ⋆-map and taking closures leads to This defines a contractive IFS on (KR) 4 , where KR is the space of non-empty compact subsets of R. This is a complete metric space when equipped with the Hausdorff metric as distance between compact sets. By the contraction principle, also known as Hutchinson's theorem in this context, the IFS has a unique solution [10, Thm. 1.1]; see [12] for related results.
One can verify that this solution is given by Since these are the windows of the original Fibonacci chain from above, we see that the disjoint point sets Λ ⋆ a and Λ ⋆ a , as well as Λ ⋆ b and Λ ⋆ b , are dense in the same sets. In particular, this shows that none of the individual point sets Λ α , with α ∈ {a, a, b, b}, is a model set. Note that the bar-swap inflation (3.1) leads to a deterministic but non-trivial partition of the a and b positions of the original Fibonacci chain into two sets each.
What is more, the images of the four types of point sets under the ⋆-map are still uniformly distributed in the corresponding windows [8,Thm. 5.3]. This implies that the FB coefficients, by the standard Weyl equidistribution argument [7,Lemma 9.4], are given by Since we know that the Bombieri-Taylor property holds for primitive inflation systems, see As the original Fibonacci system is a factor of its twisted bar-swap extension, via identifying a with a and b with b, it is clear that the dynamical spectrum of our bar-swap extension is of mixed type, which will be reflected for the generic choice of the weights u α in the diffraction measure as well. Since the extension is (measure-theoretically) 2:1, as follows from ordinary Fibonacci being an almost everywhere 2:1 factor, we know that the continuous part of the spectrum must be of pure type [5], which turns out to be singular continuous in this case, by an application of the Lyapunov exponent criterion for the absence of absolutely continuous spectral components [6,17]. The result can now be stated as follows. Let us now turn to planar systems, which will permit further variations of geometric origin.

Second variation -a planar direct product tiling
Here, we start with the obvious in taking the direct product of two Fibonacci chains, thus creating a simple tiling of the plane with four prototiles, as shown in the left panel of Figure 1. This tiling can be generated by the inflation rule where P is again a projector of rank 1. Any tiling of the hull can be viewed as a Delone set by taking the lower left corners of the tiles as control (or marker) points. In the positive quadrant, the control point sets of any member of the above 2-cycle satisfy the self-similarity relations The corresponding displacement matrix is with which the relations (4.3) simply become Λ i = j τ Λ j + T ij . By applying the ⋆-map and taking closures, this turns into where W j = Λ ⋆ j . Due to taking closures, the unions need no longer be disjoint. Eq. (4.4) defines a contractive IFS on (KR 2 ) 4 , equipped with the Hausdorff metric topology. It is easy to verify via an explicit computation that the unique solution -as expected -is given by The W j can now be interpreted as the windows (or rather the closure of the windows) of the description of the fixed points as particular projection sets, namely as regular model sets.
Clearly, for the direct product structure considered here, they are given as products of the corresponding windows for the two letters in the projection description of the Fibonacci chain; see the left panel of Figure 4 below for an illustration. Consequently, the FB coefficients or amplitudes of the defining Delone sets (for either choice from the 2-cycle) have product form and are given by in terms of the one-dimensional amplitudes from Eq. (2.1), with k = (k 1 , k 2 ) ∈ R 2 . The internal Fourier matrix is B(y) = } δ T ⋆ (y), which explicitly reads This defines the internal cocycle via B (n+1) (y) = B(y)B (n) (σy) for n ∈ N, with B (1) = B and B (n) (0) = M n . In complete analogy to the one-dimensional case, exists, with C(y) = |c(y) u| and |c(0) = |v , so that C(0) = P . For any k ∈ L ⊛ ×L ⊛ , which is the Fourier module (and the dynamical spectrum) of our direct product system, the FB coefficients are related to C by while they vanish everywhere else. As a consequence, the amplitudes can be computed, to arbitrary precision, from an approximation of C(y) via the internal cocycle. To check that this results in the same expressions for the amplitudes as the direct formula from the simple rectangular windows is left to the interested reader. The result can be summarised as follows.
Theorem 4.1. The Fibonacci direct product TDS (Y 2 , R 2 ), in its geometric realisation as described above, has pure point spectrum, both in the diffraction and in the dynamical sense. The Fourier module is L ⊛ ×L ⊛ and agrees with the dynamical spectrum of (Y 2 , R 2 ) in additive notation.
If we choose weights u 0 , . . . , u 3 for the four types of points, the diffraction measure is of the form γ = k∈L ⊛ ×L ⊛ I(k) δ k with intensity I(k) = Figure 2. Illustration of the diffraction measure γ of the Fibonacci direct product system with u 0 = u 1 = u 2 = u 3 = 1. The individual Dirac measures of γ for k ∈ [−5, 5] 2 ∩ L ⊛ × L ⊛ are represented by solid discs centred at the location of the peak, with an area that is proportional to the intensity. In later diffraction images, we shall only show the central part (grey square) with k ∈ [−2, 2] 2 for the purpose of convenient exposition.

Third variation -rearranging the direct product
Motivated by the DPVs from [13,14], we now modify the inflation rule by changing the image of the large square as follows, Then, the modified self-similarity relations in the positive quadrant read Here, the ⋆-map turns them into the IFS It turns out that the closed sets W ′ i that satisfy this IFS (5.2) are given by Indeed, inserting this into the IFS (5.2) reproduces the original square Fibonacci IFS (4.4), except for the equation for W 1 which becomes with different shifts. However, the compact sets W i of Eq. (4.5) also satisfy this equation, which corresponds to arranging the rescaled images of W 2 and W 3 in the opposite order (note that σ < 0); see the top left panel of Figure 5 for an illustration. This relation for the windows translates into a simple relation for the FB coefficients (or amplitudes). Indeed, for i ∈ {0, 1, 2, 3}, one finds where the second step follows from a simple change of variable calculation. As can be checked explicitly, S T maps the Fourier module L ⊛ ×L ⊛ onto itself. As a consequence of this analysis, this DPV also leads to a regular model set and thus to a system with pure point spectrum, both in the diffraction and in the dynamical sense.
Let us now modify the inflation rule by changing the image of the large square as follows, Iterating this inflation rule produces the tiling shown in the central panel of Figure 1. Again, the patch consisting of four large squares is legal (as can be seen in Figure 1), and it produces a fixed point under the sixth power of the inflation rule (5.4). One can determine the corresponding window IFS in complete analogy to above, but now finds windows with fractal boundaries; see the top left panel of Figure 7 below for an illustration. In fact, these windows are Rauzy fractals [19,21], for which one can show that their areas are well defined and agree with the areas of the windows from Eq. (4.5). Establishing this requires some explicit estimates of the covering regions on the basis of the contractive IFS, the details of which we omit here. Consequently, also this DPV results in a regular model set. Proposition 5.1. The dynamical system defined by the DPV rule of Eq. (5.4) has pure point spectrum, both in the diffraction and the dynamical sense. In particular, the dynamical spectrum agrees with that of the Fibonacci direct product system from Theorem 4.1, while the eigenfunctions, and thus also the diffraction measures, differ.
It is now a natural question how different the systems are which are produced by geometric variations of the Fibonacci direct product. To pursue this systematically, we consider all possible rearrangements of the stone inflation of each tile type. Clearly, there is nothing to rearrange in the inflation of a type-0 tile; there are two rearrangements each of the tiles of type 1 and 2, and 12 possibilities of the type-3 tile. These are shown in Figure 3, along with the labelling scheme used.
Note that we only consider rearrangements on the level of the tiles, and then always use the lower left corner of each prototile as marker or control point. Other choices are MLD with one of these, and do not change the spectral type. Altogether, we thus obtain 48 distinct inflation rules, all with the substitution matrix from Eq. (4.2). We parameterise the cases by triples (i 1 , i 2 , i 3 ) with i 1 , i 2 ∈ {0, 1} and i 3 ∈ {0, 1, . . . , 11}. In particular, (0, 0, 0) is the inflation from Eq. (4.1), while (0, 0, 1) is the one from Eq. (5.1) and (0, 0, 6) that of Eq. (5.4).
For each DPV, one can derive the window IFS in the same way as explained above. For precisely four choices of the parameters, namely (0, 0, 0), (0, 1,9), (1,0,3) and (1,1,6), one obtains the windows of Section 4 or a global translate thereof; see Figure 4. Observe that the inflation rules (0, 1,9) and (1, 0, 3) emerge from (0, 0, 0) by a reflection in the horizontal and vertical axis, respectively, and (1, 1, 6) by applying both reflections. Since the original tiling hull is reflection symmetric, these four rules define the same hull. As a consequence of our convention to always choose the lower left corner as control point, which does not preserve the reflection symmetry, it turns out (after some calculations) that the resulting windows are related by global shifts. This in particular implies that these four DPVs share the same diffraction, namely the one illustrated in Figure 2.
Beyond these, there are another 24 cases where the windows are parallelograms. They emerge from the original windows by shear transformations and shifts, as in the example (5.1) discussed previously. All vertices of the parallelograms are points in Z[τ ] 2 with simple coordinates. These 24 cases are illustrated in Figure 5. The proof for each case consists in determining the window IFS followed by a verification that the corresponding set of quadrangles satisfies the IFS. We note that window systems with different slopes cannot produce  Snapshots of the diffraction measures for the 12 cases with horizontal window boundaries are shown in Figure 6. Let us make some brief comments on their structure, and that of the other 12 cases. Directions in the Fourier module are mapped to directions in internal space under the ⋆-map, which is totally discontinuous. The directions of the 'white streets' (horizontal or vertical) are orthogonal to the vertical or horizontal edges of the windows. The second direction of the window boundaries shows up as a second direction in the intensity distributions, but in a less obvious way due to the complicated behaviour of the ⋆-map.
The remaining 20 DPVs lead to windows of Rauzy fractal type. Rauzy fractals are compact sets that are topologically regular and perfect. Moreover, they have positive measure and a boundary of fractal (or Hausdorff) dimension less than that of ambient space, in this case d H < 2; see [21] for a summary of the general theory. They come in three topological types, which we will call 'castle', 'cross' and 'island' to capture their geometric appearance.
There are four DPVs with windows of type 'castle' (top row of Figure 7). Two of them (the outer ones) possess a reflection symmetry with respect to the diagonal, in the sense that the windows for the control points of the squares are symmetric, while the other two are interchanged. The remaining two DPV window systems form a 2-cycle under this reflection, with the corresponding exchange of the control points of the rectangular tiles. These relations between the DPVs can directly be extracted from the corresponding inflation rules as well.
The window systems of type 'cross' and 'island' are also shown in Figure 7. There are eight in each case, whose interrelations under reflections can be studied on the level of the inflation rules. Since the windows encode the control points, which are mapped to new positions under reflections, such reflections result in more complicated rearrangements of windows and parts of windows. We leave further details, in particular a complete orbit analysis under reflections and rotations as well as the ensuing classification into MLD classes, to the interested reader.
By an inspection analogous to the one that led to Proposition 5.1, it is clear that all these examples lead to regular model sets, some of which have windows with fractal boundaries. Each individual case comes with its characteristic internal Fourier matrix, which gives rise to the matrix function C(y) = |c(y) u|, where u| and |c(0) = |v are always the same, as is the Fourier module. Since A i (k) = τ 2 5 c i (k ⋆ ), we can compute the FB amplitudes, and hence the diffraction, for all our examples via the internal cocycle. Due to the exponentially fast convergence of the cocycle product, the displayed results are free from numerical artefacts. Three examples are shown in Figure 8.
Theorem 5.2. The 48 inflation TDSs that emerge from the above DPVs all have pure point dynamical spectrum, namely L ⊛ × L ⊛ . These systems are thus metrically isomorphic by the Halmos-von Neumann theorem. Each individual tiling, via the control points, leads to a Dirac comb with pure point diffraction measure.
It remains an interesting problem to understand the topological differences between the various DPV classes, which manifest themselves in different intensity distributions of the diffraction measures, or, equivalently, in different eigenfunctions for the (commuting) Koopman operators. In particular, while the 24 examples of Figure 5 group into 12 different MLD classes, it is plausible that they all are topologically conjugate to one another, and to the Fibonacci direct product TDS. In contrast, it is clear that the Rauzy fractal windows correspond to new classes under topological conjugacy, the details of which remain to be analysed.
Clearly, the above analysis is of exemplary type and can be applied to other Pisot substitution systems as well. In particular, while we chose planar examples for ease of presentability, the method works in higher dimensions as well. In view of the increasingly frequent appearance of Rauzy fractals as windows, the internal cocycle should prove to be a valuable tool to gain further insight.