Emergence of quasiperiodic Bloch wave functions in quasicrystals

We study the emergence of quasiperiodic Bloch wave functions in quasicrystals, employing the one-dimensional Fibonacci model as a test case. We find that despite the fact that Bloch functions are not eigenfunctions themselves, superpositions of relatively small numbers of nearly degenerate eigenfunctions give rise to extended quasiperiodic Bloch functions. These functions possess the structure of earlier ancestors of the underlying Fibonacci potential, and it is often possible to obtain different ancestors as different superpositions around the same energy. There exists an effective crystal momentum that characterizes these ancestors, which is determined by the mean energy of the superimposed eigenfunctions, giving rise to an effective dispersion curve. We also find that quasiperiodic Bloch functions do emerge as eigenfunctions when weak disorder is introduced into the otherwise perfect quasiperiodic potential. These theoretical results may explain a number of experimental observations, and may have practical consequences on emerging theories of band topology and correlated electrons in quasicrystals.


I. INTRODUCTION AND MOTIVATION
Quasicrystals are quasiperiodic crystals that are strictly aperiodic [1]. Their discovery [2] put an end to the age-old paradigm that long-range order-in the positions of atoms in a material-is synonymous with periodicity, ushering in a Kuhnian scientific revolution [3]. In the past four decades we have experienced the exciting paradigm-rebuilding phase of this revolution [4], in which old notions and ideas are being reexamined, and carefully modified and adapted to the age of quasicrystals [5]. Of particular relevance here is the fundamental notion that order, or lack thereof, is what determines the nature of single-particle quantum excitations in condensed matter.
Bloch [6] showed, already in 1928, that electrons in periodic potentials form extended eigenfunctions, almost as if they were in free space, with absolutely continuous energy spectra. Anderson [7] established three decades later that electrons in disordered potentials possess exponentially localized eigenfunctions, with pure point spectra. Surprisingly, only in dimensions greater than two is there a delocalization transition as a function of energy [8], whereby beyond a critical value-the so-called mobility edge-electrons have sufficient energy to overcome the disorder for their eigenfunctions to become extended. It was only natural for the scientific community, which at the time was still unaware of aperiodically ordered matter, to associate the existence of extended or localized eigenfunctions with having ordered or disordered potentials, respectively. The discovery of quasicrystals called for a reexamination of this notion.
Quasicrystals possess long-range order. It is therefore natural to expect their single-particle eigenfunctions to be extended and obey a generalized form of Bloch's theorem for quasiperiodic potentials. Yet, four decades of * Corresponding author: ronlif@tau.ac.il research on single-particle excitations in quasicrystals  have shown that this is not the case. Depending on the nature of the model used, single-particle eigenfunctions in quasicrystals can take almost any possible form: they can be extended as in periodic potentials, localized as in disordered potentials, or critical, exhibiting powerlaw or algebraic decay, which is not seen in the former cases.
As it seems, the Bloch theorem, which is the cornerstone of the theory of electrons in periodic crystals, does not generally apply to aperiodic crystals. Rather than proposing alternative forms for the electronic eigenfunctions in quasicrystals [30,31], or suggesting a generalized formulation for the Bloch theorem [32], we examine here whether Bloch wave functions may still appear naturally in quasicrystals even when Bloch's theorem fails to generate them as eigenfunctions. This may have immediate practical implications, as the study of correlated electrons in quasicrystals [33,34] is regaining interest [35][36][37], amid the experimental observation of Dirac electrons in dodecagonal bilayer graphene quasicrystals [38,39], and the reformulation of the early spectral theories of electrons, phonons, and photons in quasicrystals [40,41] using modern topological classifications [42][43][44][45][46]. It may also have experimental significance with the increasing ability, at least in metamaterials, to observe and manipulate the actual wave functions [47] and their topological nature [48][49][50], while inducing topological defects [51], phason excitations [52], and static disorder [53] in the otherwise perfectly ordered quasicrystalline potential.
To be concrete, let us consider a quasiperiodic crystalline potential U (r), which can be decomposed into countably many Fourier modes, where the reciprocal lattice L, consisting of the closure under addition of all wave vectors k with nonzero Fourier arXiv:2108.11750v2 [cond-mat.mtrl-sci] 26 Mar 2022 coefficientsŨ (k), is a finitely generated Z-module. If the rank, or smallest number D of wave vectors required to generate L over the integers, is equal to the spatial dimension d, then U (r) is the potential of a periodic crystal. In other words, if D = d, L is reciprocal to a Bravais lattice T of real-space translations leaving the crystalline potential invariant, in the sense that U (r + R) = U (r) for every R ∈ T . If D > d, then U (r) is the potential of a quasicrystal, there exists no lattice of spatial translations leaving the potential invariant, and the reciprocal lattice L becomes dense and no longer possesses a Brillouin zone. According to the Bloch theorem (see, for example, chapter 8 of Ashcroft and Mermin [54]), single-electron eigenfunctions, satisfying the Schrödinger equation for such a potential, can be expressed as where the quantum number q labeling them is known as the crystal momentum. Consequently, the Fourier transforms of the electron density |ψ q (r)| 2 and the potential U (r) are supported on the same reciprocal lattice L, while the Fourier transform (2) of the eigenfunction itself is supported on the shifted lattice L + q = {k + q | k ∈ L}. By comparing the Fourier transforms of |ψ q (r)| and ψ q (r), one can readily determine the crystal momentum q by the presence of a uniform shift between the two spectra. Repeating this process for all eigenfunctions yields a dispersion relation E(q), which corresponds to the crystalline band structure. This holds true if the crystal is periodic. It generally fails if the crystal is aperiodic, possessing a dense reciprocal lattice, owing to the fact that the infinite sum over scattered waves in Eq.
(2) generically does not converge, even though the corresponding sum in Eq. (1) does. Nevertheless, it is interesting to explore whether extended quasiperiodic Bloch functions may somehow emerge in realistic physical situations-either as naturally occuring superpositions of eigenfunctions, or even as true eigenfunctions of slightly modified quasiperiodic structures, for example, by introducing static disorder or external fields. If so, what would be the nature of these quasiperiodic Bloch functions? How would they be related to the structure of the underlying quasicrystalline potential? Could they be characterized by similar quantum numbers q in reciprocal space as their periodic analogs? If so, would there exist an energy-momentum dispersion curve, or effective band structure, that could be associated with these Bloch functions? Positive answers to these questions may explain a number of empirical observations that seem to indicate that energymomentum dispersion curves do exist [55,56], and that slight disorder may increase the conductivity [57] and the spatial extent of wave functions [53]. The latter effect has been demonstrated recently using renormalization-group calculations [58].

II. CHOICE OF MODEL
We would like to choose, among the commonly used models for studying electrons in quasicrystals, one that is as simple as possible, yet sufficiently generic to explore the questions raised above. We therefore concentrate on one dimension-although some initial calculations of ours in two dimensions can be found elsewhere [59,60]-and limit ourselves to quasicrystals of rank D = 2, that is with just a single pair of incommensurate fundamental spatial harmonics. With these restrictions there are still a few different types of models to choose from, some more physically motivated than others. These can roughly be categorized as follows (with more details in the references provided):

A. Continuous time-independent Schrödinger equations
These are ordinary differential equations of the form with continuous quasiperiodic potentials U (x) as in Eq.
(1), supported on a reciprocal lattice L of rank D = 2 [9][10][11][12]. To simplify things, one often limits the potential to its two fundamental harmonics by taking, for example, with an irrational τ . For weak coupling λ 1, or weak quasiperiodicity 1, one can treat the equation as a small perturbation with respect to free electrons, or to Bloch electrons in a periodic crystal, respectively. In these limits the spectrum remains purely absolutely continuous, exhibiting a well-defined band structure E(q), with a hierarchy of gaps that open at q = nπ + mπτ (n, m ∈ Z), as a result of the hybridization that occurs whenever the degenarate ±q free-electron states differ by a wave vector k ∈ L. Formally, this happens because of small divisors that appear in the perturbation expansion (see, for example, chapter 9 of Ashcroft and Mermin [54]). In these limits the eigenfunctions are quasiperiodic and satisfy Bloch's theorem (2).
Owing to the dense nature of L, for sufficiently large λ and , and depending on the Diophantine properties of the irrational ratio τ , the formation of gaps typically destroys the band structure at the bottom of the spectrum. Consequently, the bottom of the spectrum becomes pure point with eigenfunctions that are exponentially localized. Typically [61], there is a critical energy E C = E C (λ, , τ ) below which the eigenfunctions are exponentially localized and above which they are extended, akin to the mobility edge that is observed in Anderson localization in dimensions greater than 2. It should be emphasized, though, that while Anderson localization arises from the lack of order, the localization here arises from the existence of order, albeit aperiodic order. It is the long-range order that is responsible for having strong Bragg peaksŨ (k) in the Fourier transform of the potential (1), that can then combine with the dense nature of L to destroy the continuous bands.

B. Discrete time-independent Schrödinger equations
These are finite difference equations of the form with a periodic analytic function f (x) whose period is incommensurate with that of the lattice m ∈ Z on which it is sampled. A standard example is to take with an irrational τ . In these finite-difference, or tight-binding, Schrödinger equations (5), known as the Aubry-André model [13], or the almost Mathieu equation [11,14,[62][63][64], one no longer observes a localization transition as a function of energy. Instead, what one typically finds is that there is a critical value λ C of the coupling constant, equal to 2 for the standard example (6), at which the whole spectrum changes its nature [64]. For almost every choice of irrational τ and real θ, if λ < λ C the bands remain intact and the whole spectrum is purely absolutely continuous with eigenfunctions that are extended. If λ > λ C the whole spectrum is pure point and the eigenfunctions are exponentially localized. Exactly at λ = λ C the spectrum is purely singular continuous [63], like a Cantor setuncountable yet of zero Lebesgue measure-neither pure point nor consisting of any bands or continuous intervals. The corresponding eigenfunctions decay algebraically.

C. Tight-binding time-independent Schrödinger equations on quasiperiodic tilings
In one dimension, quasiperiodic tilings reduce to quasiperiodic sequences of finitely many different letters. These letters may then represent either the sequence of distances between atoms on the line, affecting the hopping amplitudes {T m } between neighbors, the sequence of different atomic species arranged on the line, affecting the onsite energies {ε m }, or both, as given by When the T m are all equal, or the ε m are all equal, Eq. (7) reduces to the diagonal or the off-diagonal tight-binding models, respectively.
The essential difference between these equations and those of the previous category is that the sequence f (m) in Eq. (5) consists of infinitely many different values, densely sampling the image of the continuous function f (x), whereas the sequences {T m } or {ε m } consist of a finite number of distinct values. A standard example is the family of Sturmian sequences [65], consisting of two letters, corresponding to whether vertical or horizontal lines are crossed when cutting through a square grid with an irrationally sloped straight line. This includes the most familiar example of the Fibonacci sequence [15][16][17][18][19][20], recently reviewed by Jagannathan [66].
The striking behavior for a large family of these models [65][66][67][68][69] is that they exhibit purely singular-continuous zero-Lebesgue-measure spectra, with critical eigenfunctions, for any nontrivial choice of their parameters. They never exhibit absolutely continuous spectra with extended Bloch eigenfunctions, as do periodic crystals, and they never exhibit pure point spectra with exponentially localized eigenfunctions, as disordered structures do. They represent what is most unique about quasicrystals in this respect, and are therefore best suited to study the main question posed here: Is it still possible for Bloch wave function to occur naturally in these aperiodically ordered models, despite the failure of Bloch's theorem in any part of their spectra and for any choice of their parameters.

D. Time-independent quasiperiodic Kronig-Penney models
The last family of models that should be mentioned in this context are quasiperiodic generalizations of the 1931 Kronig-Penney model [70]. These models [22][23][24][25][26][27][28][29] are in some sense intermediate between the continuous models of Sec. II A and the discrete models of Secs. II B and II C, as they replace the continuous potential U (x) in the Schrödinger equation (3) by a discrete sum of delta functions with variable weights W m , These quasiperiodic Kronig-Penney models typically exhibit the same features as the discrete tight-binding versions, usually with some added complexities that are lost in the tight-binding approximation. An interesting example is the appearance of countably many extended Bloch eigenfunctions, whose energies are embed-ded within the singular continuous spectrum of an uncountable Cantor set of critically decaying eigenfunctions [28,29].
E. Our choice: The off-diagonal tight-binding Fibonacci model As argued in Sec. II C above, we choose to consider the family of tight-binding models of Eq. (7), where the Bloch theorem fails completely, and the eigenfunctions are all critically decaying for any nontrivial choice of parameters. We employ the well-known Fibonacci sequence, taking the off-diagonal version of the model for simplicity, Here, the onsite energies are all zero, and the hopping amplitudes Accordingly, the length of the finite sequence S n is the Fibonacci number F n , where F n = F n−1 + F n−2 and F 1 = F 2 = 1. This is the same as starting with the sequence S 1 = {S} and iteratively applying the substitution rules S → L and L → LS. Moreover, the finite sequences S n yield the optimal choices of unit cells for periodic approximants of the infinite Fibonacci quasicrystal. A finite sequence S m is said to be an ancestor of the sequence S n if m < n.
The approximant tight-binding Hamiltonian, corresponding to Eq. (9), is thus given by a F n × F n matrix that is tridiagonal up to boundary terms. A direct diagonalization of the matrix, for a given value of T and a particular choice of boundary hopping terms, is then used to obtain a discrete set of F n eigenvalues along with their corresponding eigenfunctions, converging to the quasicrystalline spectrum in the limit of n → ∞. A typical spectrum is shown in Fig. 1 for an approximant with N = F 16 = 987 sites, and T = 6. Typical algebraically decaying eigenfunctions are shown in the bottom inset of Fig. 1, as well as in blue in the top three panels of Fig. 2, for an approximant with N = F 15 = 610 sites, and T = 6. Indeed, for any T = 1 the eigenfunctions all tend in the infinite limit to critical wave functions that decay algebraically [15][16][17], rather than to extended Bloch functions. FIG. 2. The top three panels show three typical spatially decaying eigenfunctions of the Fibonacci quasicrystal, with nearly degenerate consecutive eigenvalues, as indicated within the panels, differing by less than 10 −5 of the total spectral width. The bottom panel shows a linear combination of these three eigenfunctions, whose peaks are spatially extended, and for which the mean energy is E = −5.656905, with an uncertainty, or standard deviation, of ∆E = 4.2 · 10 −5 . Here T = 6 and N = F15 = 610.

III. LINEAR COMBINATIONS OF NEARLY-DEGENERATE EIGENFUNCTIONS
We wish to see whether extended Bloch wave functions may still appear naturally in these tight-binding models, albeit not as eigenfunctions. As a first step, and following early ideas of Even-Dar Mandel and Lifshitz [59,60,71], we consider the spatial extent one can achieve by taking linear combinations of nearly degenerate eigenfunctions. This idea is motivated by the fact that many of the simple structural or tiling models of quasicrystals are linearly repetitive. This means that any finite patch of radius r is repeated in the tiling at a distance that scales linearly with r. Thus, algebraically decaying eigenfunctions that are very close in energy, which are likely to originate from large patches that are structurally similar, will have their peaks located at distant positions in the quasicrystal. This may allow a relatively small number of nearly degenerate eigenfunctions to span the whole quasicrystal. We note that Niu and Nori [17] already saw indications that this might be the case by considering incoherent sums |ψ n | 2 of nearly degenerate eigenfunctions, but not as coherent wave functions c n ψ n that may describe an actual particle with a small uncertainly in its energy.
The three eigenfunctions shown in blue in the top three panels of Fig. 2 correspond to three consecutive eigenenergies in the spectrum of a N = F 15 = 610 approximant, and together may indeed form a coherent superposition that spans the whole approximant, as shown in red in the bottom panel. The perceptive reader may notice that the peaks in the extended function are separated by long (L) and short (S) intervals, arranged according to the N = F 6 = 8 Fibonacci approximant LSLLSLSL. This reflects the spatial distribution of similar local sequences, in the F 15 approximant, inherited from its ancestor of 9 generations earlier.
To pursue this intuition, we numerically optimize the spatial extent of linear combinations of nearly degener- ate eigenfunctions, whose energies lie within tiny windows W around every eigenenergy E n in the spectrum of a given approximant, as demonstrated schematically in the top inset of Fig. 1. It is most often the case that several consecutive energies yield the same result owing to the near-flatness of the spectrum. We use the difference-map algorithm of Elser et al. [72], with the aim of generating a wave function that satisfies two constraints simultaneously: (I) it should be as extended as possible, ideally having an equal amplitude on all sites; and (II) it should be spanned by the assigned set of nearly degenerate eigenfunctions. The two constraints are applied iteratively-the first by setting |ψ(m)| = 1/ √ N on all sites, while keeping the phases, and the second by projecting the wave function into the subspace spanned by the given eigenfunctions-always finishing with the latter projection [59,60]. We note that more traditional constraint-solving approaches, like least-squares algorithms, yield similar results, albeit with longer computation times. To assess the results and the progress of the calculation, we use the Normalized Participation Ratio, as a measure of extent. By construction, the NPR is equal to 1 for a uniformly extended wave function, and decays to 0 as it becomes localized. Surprisingly, the wave functions that emerge from this optimization all happen to be Bloch-like Fibonacci wave functions, even though they are merely optimized for extent. They are Fibonacci wave functions in the sense that they have the structure of an earlier Fibonacci ancestor of the underlying potential, as anticipated in the bottom panel of Fig. 2. They are Bloch-like wave functions in the sense that their Fourier spectra, as shown in Fig. 3, are carried by the same set of wave-vectors as the underlying quasiperiodic potential-set by the Fibonacci sequence of hopping amplitudes T m -shifted by the crystal momentum q.
Although we are limited computationally up to approximants of size N = F 18 = 2584, it seems that this behavior persists with increasing N . Specifically, we find that while the NPRs of the individual eigenfunctions decrease with N , the NPRs of the most extended linear combinations of eigenfunctions remain roughly constant with increasing N , suggesting that the Bloch-like Fibonacci character of the extended linear combinations may persist in the thermodynamic limit.
Moreover, it is possible to extract the crystal momentum q numerically for each energy, by comparing the Fourier transform of the potential to that of the extended wave function. Specifically, we calculate q by the shift in the main, or strongest, peaks in the Fourier spectra of the wave function and the potential. We note that the accuracy in which we are able to determine q improves with increasing approximant size. However, if the most extended wave function is a very early ancestor, as is the case in spectral regions with a low density of states, the resemblance between the Fourier spectra is weaker, and the determination of q becomes less accurate.
The resulting effective dispersion relation E(q), displayed as red diamonds in Fig. 4, is very similar to the eigenenergy spectrum in Fig. 1, suggesting that q might be proportional to n. This is reminiscent of the situation for periodic crystals, where q = 2πn/N , although some overlap exists between the apparent minibands, spoiling the perfect monotonic dependence of E on q. We are unable to determine whether this overlap is real, or a result of not having obtained the optimal linear combination, leading to an inaccurate mapping of q to E.
Looking more quantitatively, we find that both the NPR and the ancestral generation, or Fibonacci number, associated with the most extended wave functions, are correlated with the spectral density of states around the mean energy of the extended wave function. As the density of states increases, a given fixed energy window W will contain more eigenfunctions, allowing one to obtain a more extended wave function, and to more closely follow the Fibonacci potential, yielding a more recent ancestor. With values of T large with respect to unity, as is the case for T = 6, shown in Fig. 1, the spectrum appears as a hierarchy of rather flat clusters of eigenvalues-each corresponding to a Fibonacci number of eigenfunctionsand the density of states is sharply peaked. In such cases a very small window W , of less than 10 −3 of the full bandwidth, already contains a fairly large cluster of 21-55 of the 987 eigenfunctions, and a highly extended wave function corresponding to a recent ancestor is obtained. As W decreases, fewer functions participate in the linear combination, the NPR of the optimized wave function gradually decreases, and the ancestral order increases, while the quantum number q changes very little, although its numerical determination becomes less accurate. This qualitative behavior remains similar for all values of T . One should only note that as T approaches unity, deviations from periodicity decrease, gaps in the spectrum gradually close, and the NPRs of the eigenfunctions themselves increase. Rather than optimizing for extent, we can employ standard distance-minimization algorithms using earlier ancestors as targets for optimization, without decreasing the number of eigenfunctions used. The result of such a construction is shown in Fig. 5, where four different ancestors are obtained as linear combinations of the same set of nearly degenerate eigenfunctions. We find that these earlier ancestors have similar values of q, which change monotonically with their mean energies. This seems to imply that earlier ancestors provide additional samples of the dispersion relation E(q) = ψ q | H |ψ q , generating a denser dispersion curve, shown as blue dots in Fig. 4. However, this curve exhibits more overlap between the minibands, with q ∝ n only within each miniband, perhaps because the determination of q becomes less accurate the earlier the ancestor.

IV. ADDITION OF WEAK DISORDER
Our analysis indicates that a small uncertainty ∆E W , in the energy of a single electron, may allow it to ex- plore sufficiently many nearly degenerate energy eigenfunctions for it to behave like a Bloch electron. Such an uncertainty may arise naturally from semiclassical dynamics in a "semiadiabatic" regime [73] that smooths out small gaps in the spectrum, leaving effective minibands E(q) between the remaining large gaps. Alternatively, as we explore here, a mixing of nearly degenerate eigenfunctions may be induced by adding weak disorder to the otherwise perfect quasicrystalline Hamiltonian (9). We model the disorder using zero-mean Gaussian random variables with standard deviation ∆T , added to the hopping amplitudes T m in Eq. (9). Similar results are obtained for δ-correlated, or white, disorder as for inverse power-law correlated disorder, so we focus here on the former. We perform an ensemble average of typically 1000 disorder realizations for each eigenfunction. The disorder is kept sufficiently weak to maintain monotonic dependence of the energy eigenvalues on disorder strength, so as not to change the identity of the eigenfunctions by approaching level crossings.
As expected, extremely weak disorder hardly affects the eigenfunctions, while sufficiently large disorder leads to Anderson localization, i.e. to exponentially decaying eigenfunctions. A little less expected, yet consistent with previous observations [53,57,58], is the result that weak disorder, on the order of the windows W used earlier or less, improves the extent of some of the eigenfunctions, while reducing the extent of others. This is demonstrated in the histograms of Fig. 6, where one can see that the sharply distributed NPR values in the perfect Hamiltonian are more evenly distributed in the ensembleaveraged disordered Hamiltonian. We note that for the disorder strength giving rise to these results, the energy spectrum itself remains nearly unchanged.
More important is the observation that weak disorder significantly increases the NPR of the most extended eigenfunctions, and that the resulting ensemble-averaged functions take the form of nearly perfect early Fibonacci ancestors, like the ones obtained above by optimizing linear combinations of eigenfunctions of the perfect Hamiltonian. However, since these are very early Fibonacci ancestors, we are unable to faithfully determine their corresponding wave vector q. The most extended eigenfunctions, obtained with added disorder of a given strength, are quite similar to each other and resemble the Fibonacci ancestor shown in Fig. 7.

V. CONCLUSION
The study of quasicrystals has taught us that the electronic properties of materials depend crucially, not only on how ordered they are, but also on how they are ordered. As we demonstrate here, the interplay between disorder and aperiodic long-range order is particularly intriguing.
Owing to the dense nature of the Fourier module L, the Bloch sum over wave vectors in Eq. (2) is not guaranteed to behave properly for quasicrystals, as it does for periodic crystals. As a consequence, Bloch's theorem generally fails and eigenfunctions are generically critical, decaying algebraically rather than being quasiperiodically extended throughout the crystal. Nevertheless, quasiperiodic Bloch wave functions do form as superpositions of relatively small numbers of nearly degenerate critical eigenfunctions, which emerge naturally in the presence of weak disorder. Contrary to its effect on periodic crystals, disorder in quasicrystals first increases the extent of the most extended eigenfunctions, transforming them from critical to extended quasiperiodic Bloch functions, before eventually giving way to Anderson's inevitable localization.
As for periodic crystals, one can associate an effective crystal momentum q with each quasiperiodic Bloch wave function that forms, leading to an effective dispersion relation E(q), which may explain the dispersion curves observed in certain experiments [55,56]. We expect this behavior to occur in other linearly repetitive quasicrystalline models, and to be even more pronounced in two and in three dimensions, where Anderson localization is less restrictive. We leave the verification of this expectation, as well as the analytical explanation of our empirical findings, and their extension to other models-perhaps considering models like the Thue-Morse sequence that are deterministic yet possess no long-range order-for future research.