Searching for metal-insulator transition in Fibonacci quasicrystals by distribution analysis of electronic wavefunction participation ratios

Since their discovery in 1984 by Shechtman [1], the electronic properties of quasicrystals (QCs) have always attracted special attention, as their quantum-mechanical wavefunctions can actually be considered intermediate in character between the extended states of a metal and the exponentially localized states of an insulating phase. Several authors analyzed the electronic properties of such systems, taking into account many-body effects in the Fibonacci chain [2] and the metal-insulator transition in two-dimensional (2D) QC from the Amman-Beenker tiling family in the framework of Random Matrix Theory (RMT) [3]. Investigations concerning the distorsion-induced metal-insulator transition in simpler tight-binding models for 1D or 2D QC tilings from the metallic-mean family [4] have not been published so far in literature. It would be far more interesting to explore the likely connection between the metallic number [5] defining the QC aperiodic structure and the onset of metal-insulator transition, not only on a fundamental science viewpoint but also with the aim of the designing new technological applications based on metamaterials [6] where this phenomenon can be controlled and eventually tuned.

The present note is focused on the analysis of the shape change in the distribution of participation ratios of the computed QC electronic wavefunctions when the geometric aperiodic structure is distorted in order to drive a phase transition in the electronic spectrum of the system. The normalized participation ratios are defined accordingly to equation (7.7) in [7] in terms of eigenvector coefficients e_i (E_ν ) corresponding to states e_ν with energy E_ν

$$\begin{eqnarray}&&{p}_{\nu }=\displaystyle \frac{{\left({\displaystyle \sum }_{i}| {e}_{i}({E}_{\nu }){| }^{2}\right)}^{2}}{N{\displaystyle \sum }_{i}| {e}_{i}({E}_{\nu }){| }^{4}}\end{eqnarray} \tag{ 1 }$$

where the sum over the orbital characters of the wavefunction α in the original reference has been dropped, since the investigation refers to a simplified model with a unique orbital character. The sums over index i appearing in the above expression run over all the N sites of the system. The (normalized) participation ratios p_ν quantify the degree of localization of an electronic wavefunction, namely p_ν ≪ 0. 5 points out for a localized wavefunction, whereas for p_ν ≳ 0. 5 it is extended. The electronic wavefunctions are numerically computed as the eigenvectors of the nearest-neighbor tight-binding hamiltonian matrix spanning the quasiperiodic lattice approximant of the system

$$\begin{eqnarray}&&{H}_{j,k}=\displaystyle \sum _{j,k}{t}_{j,k}\left({c}_{j}^{\dagger }{c}_{k}+{c}_{j}{c}_{k}^{\dagger }\right)+\displaystyle \sum _{k}{\epsilon }_{k}{c}_{k}^{\dagger }{c}_{k}\end{eqnarray} \tag{ 2 }$$

where ${c}_{j}^{\dagger },{c}_{k}$ are the second quantization electron creation and annihilation operators at sites j and k, respectively. ${t}_{j,k}$ gives the hopping matrix elements between electron states associated with the corresponding pair of sites. In a sparse-matrix representation of the hamiltonian graph, the off-diagonal elements given by ${t}_{j,k}$ are non-zero for vertices connected by an edge, otherwise ${t}_{j,k}=0$. The actual values of ${t}_{j,k}$ depend on the length of the edge connecting two neighboring vertices. The on-diagonal parameter _k correspond to the on-site energies, but usually the simplification _k  = 0 is adopted, as chosen here, so that no vertex is energetically favoured. This is the simplest model able to capture the essential electronic features of the system, since the aperiodicity is directly embedded in the graph modeling the QC patch and in the modulation of the hopping probabilities according to the distance between neighboring vertices.

For an ideal aperiodic crystal the matrix to be diagonalized would be infinite; in practice, as chosen in this case, one usually approximate it with a finite QC patch with both a statistically significant and computationally manageable number of sites (i.e. ∼10⁴) to extrapolate results valid for a large number N of sites. The simplest 1D quasicrystal models for patch construction are based on the Fibonacci quasiperiodic word sequence (see figure 1), which is obtained starting from a two-letter alphabet (L S) and a substitution rule σ(L) = LS, σ (S) = L that can be iterated n times to give a reasonably large QC system [6]. The ratio of the frequencies of the two letters S and L in the limit of n tending to infinity is the golden number $\tau =(1+\sqrt{5})/2$, which is the first in the series of the metallic means family [5]. Therefore, by considering equation (2), the hopping matrix elements ${t}_{j,k}$ can either take the value assigned for L or S.

Zoom In Zoom Out Reset image size

Figure 2 shows a finite patch of a 2D Fibonacci QC with a square grid obtained by performing the tensorial product of two identical 1D Fibonacci quasiperiodic sequences of different spacings L (long) and S (short) on the grid axes.

Zoom In Zoom Out Reset image size

To address the question of change in the electronic structure at the origin of the metal-insulator transition, the effects of the change of the L/S ratio on the Hamiltonian were investigated in a series of computations performed by setting L = 1 and varying S in the range 0 < S < 1 (at S = L = 1 the system is a periodic crystal) and then subsequently by letting S = 1 and 1 < L ≤ t₁ with t₁ ≫ τ, where the measure units are to be considered adimensional, because of the theoretical nature of the considered system. Histograms of the normalized participation ratios p_ν of the full set of the numerically computed tight-binding eigenvectors were collected for each L, S values in the above mentioned ranges with bin size set to 0.1.

It can be verified that the histogram distribution of normalized participation ratios obtained for L = 1, S = τ⁻¹ ≈ 0.618 matches the one obtained from S = 1, L = τ ≈ 1.618 , as a consequence of the symmetry imposed by the Fibonacci substitution rule (figure 3). The change in the p_ν distribution as a function of the L/S ratio can be more easily appreciated in the dynamical histogram plot reported in figure 4(a). When $L=1,S\ll {\tau }^{-1}$, most p_ν values concentrate on the left side of the histogram with 0.1 ≲ p_ν  ≲ 0.3, meaning that all the wavefunctions are localized [7].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By increasing S, a change in the histogram shape occurs around L = 1, S = τ⁻¹, with the largest fraction of p_ν values concentrating now around 0.5, i.e. wavefunctions are now stepping into the delocalized or extended regime. For $L=1,\ {\tau }^{-1}\ll S\lt 1$, the system is completely in the metallic regime, since all participation ratios now accumulate on the rightmost part of the histogram with p_ν  ≃ 0.7. This situation is mirrored in the range $S=1,\ L\gg \tau$, where the p_ν bins start to move back from right to left, with a change occuring again at L = τ. From that point on (leftwards), the system enters again the localized regime with $S=1,\ L\gg \tau$. This fact was reproducibly verified for several increasing n values in the substitution sequence. Here only the results for n = 20 system were reported, which correspond to a matrix size of 17 712, the largest which could be straightforwardly diagonalized on an ordinary laptop computer and a working Python software installation for numerical computing. A similar result is found for the 2D Fibonacci QC patch, as reported in figure 4(b) for n = 9 (corresponding to a matrix size equal to 8100).

According to the definition of unnormalized participation ratios ${P}_{\nu }={{Np}}_{\nu }$, which quantifies the number of lattice sites over which the wavefunction is extended (i.e. $1\lesssim {P}_{\nu }\lesssim N$), the full set of P_ν values for the 1D Fibonacci QC with n = 20 was tentatively fitted against the Polya-Aeppli (compound geometric Poisson) probability distribution with shape parameters r and q ¹ .

The plots reported in figure 5 show the r, q values resulting from distribution fitting with the Wolfram Mathematica^® function ${\mathsf{FindDistributionParameters}}$ using the builtin ${\mathsf{MethodOfMoments}}$ option for different L, S settings corresponding to the localized regime (top left, bottom right panels), the critical ratios at τ⁻¹, τ (top and bottom center panels) and the extended regime (top right, bottom left panels).

Zoom In Zoom Out Reset image size

Figures 6(a), (b) summarize the corresponding trends of r and q parameters reported on a logarithmic scale as a function of the L/S ratio in the 1D and 2D QCs, respectively. Both plots show the divergence of r, q values at the quasiperiodic-to-periodic transition at L = S = 1 marked by q = − 1. The q parameter is very close to 1 and has a very weak dependence on L/S ratio for the 1D system. The values of the size parameters r and q at τ⁻¹ and τ are highlighted with dashed vertical lines. The turning point close to r = 20 is also marked by a dashed horizontal line as a guide to the eye for both systems. By numerically interpolating the r values as a function of the corresponding L/S ratios, a flex in such functional dependence can be found in the proximity of the critical ratios (τ⁻¹, τ), which can be interpreted as the metal-insulator transition point for the 1D system. For the 2D QC a small dip is found at τ². In the latter case the computed r, q values refer to a global average over the whole set of participation ratios, without discriminating the P_ν values according to the QC orthogonal axes x and y.

Zoom In Zoom Out Reset image size

A weakly oscillating behaviour of the r values as a function of the QC size n can be seen in figure 7, which shows that the electronic properties of the 1D QC system are symmetry dependent on n parity. However, the numerical differences in the values of r according to the n parity become smeared out after n = 20 where both r trends pair asymptotically while slowly continuing to decrease.

Zoom In Zoom Out Reset image size

If more computing resources were available in terms of memory and time, it would be possible to obtain participation ratio datasets for larger n values, not only for the 1D system but also for the 2D Fibonacci QC, in order to address the scaling of the distribution shape parameter r with the system size n. In particular, it would be helpful to ascertain the existence of a limit value for the fitting parameter r, as reported in figure 7. Certainly, it would be far more interesting to provide a more in-depth physical interpretation for the shape parameters for other structures spanned by other metallic numbers beside the Fibonacci golden ratio. More generally, on a theoretical sideview, it would be compelling to investigate the relationship between the distribution of participation ratios in presence of disorder and the distributions of the energy level spacings known from RMT, according to the recently found linear algebra identity between eigenvalues and eigenvector modules [8].

In summary, the analysis of the participation ratios from single-particle eigenfuctions for low-dimensional Fibonacci quasicrystals suggests a metal-insulator transition related to a Peierls distorsion in a quasiperiodic tight-binding model with nearest-neighbour staggered hopping amplitudes (half-filled Su-Schrieffer-Heeger model) which likely occurs at a lattice distortion ratio corresponding to the golden number or its inverse, according to the chosen hamiltonian parametrization.

The author gratefully acknowledges Istituto d'Istruzione Superiore Cattaneo-Mattei PDIS003009 in Monselice-Conselve (Padova, Italy) for the supportive environment.

The data that support the findings of this study are available upon reasonable request from the authors.

The Python code for reproducing the results is available upon request from the author.